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Chapter 1 



Introduction 



1.1 Motivation and Problem Definition 

Traffic networks - consisting of highways, streets, and other kinds of roadways - 
provide convenient and economical conveyance of passengers and goods. The ba- 
sic activity in transportation is a trip, defined by its origin/destination, departure 
time/arrival time and travel route. A myriad of trips interact on the network to 
produce an intricate pattern of traffic flows. Since traffic conditions in many major 
metropolitan areas are becoming increasingly congested, affecting the operational ef- 
ficiency of whole networks as well as the travel cost of each trip, traffic flow models are 
becoming more important in traffic engineering and the transportation policy making 
process. For example, well-developed traffic models are used in developing advanced 
ramp metering methods as well as in determining dynamic traffic assignment (JWL 
& Zhang, 2000a). 

There have been two approaches in mathematical modeling of traffic flow. One 
approach, from a microscopic view, studies individual movements of vehicles and 
interactions between vehicle pairs. This approach considers driving behavior and 
vehicle pair dynamics. But the size of the problem in a microscopic model becomes 
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mathematically intractable when a considerable volume of traffic flow is considered. 
One example of the microscopic approach is the GM family of car-following models 
developed in the 1960's (e.g., Gazis et al., 1961). The other approach studies the 
macroscopic features of traffic flows such as flow rate q, traffic density p and travel 
speed v. The basic relationship between the three variables is: q = pv. Macroscopic 
models are more suitable for modeling traffic flow in complex networks since less 
supporting data and computation are needed. 

In this thesis macroscopic traffic flow models are studied both theoretically and 
numerically. Traffic flows are classified according to traffic conditions, roadway con- 
ditions and traffic network structure. Traffic flows are in equilibrium when the travel 
speed of these flows is uniquely determined as a function of traffic density, otherwise 
they are in non-equilibrium. Traffic flows are considered inhomogeneous when the 
roadway has different parameters at different locations. Link flows are flows on road 
links, and network flows are traffic flows on networks of roadways. Network flows dif- 
fer from link flows in that vehicles in the former have different characteristics which 
affect traffic dynamics, such as the origins or destinations. 

Different types of traffic flow are described by different models. For equilibrium 
link flow, the celebrated LWR model was developed by Lighthill and Whitham (1955) 
and Richards (1956). The LWR model has been solved for the homogeneous roadway 
rigorously. There have been empirical solutions to the inhomogeneous LWR model. In 
this work a rigorous procedure to solve the inhomogeneous LWR model is developed. 
(JWL & Zhang, 2000b). The LWR model is a first-order model in the sense of PDE 
system order. In this thesis we also discuss the PW model (Payne, 1971; Whitham, 
1974) and Zhang's model (1998, 1999a) for non-equilibrium link flow. Finally we 
introduce a multi-commodity model when traffic flow is disaggregated by origins, 
destinations or departure times. 

All the models we consider are based on conservation of traffic flow, i.e., the 
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increment of vehicles in a section is equal to the difference between upstream in- 
flux and downstream out-flux in unit time. Mathematically, every model except the 
multi-commodity model, which is a discrete model, can be written as a continuous 
system of hyperbolic conservation laws. Solutions to the Riemann problem for all the 
continuous models are studied analytically, and all the models including the multi- 
commodity model are solved numerically with Godunov type of methods. 

1.2 Background and Research Overview 

1.2.1 Traffic Flow Models 

In what follows the term "traffic flow model" means a macroscopic traffic flow model. 
Many continuum traffic flow models can be described by a system of hyperbolic PDEs. 
The first of these models was the LWR model. This model relies on the assumption 
that there exists an equilibrium speed-density relationship v = f*(p). Like other 
dynamic continuum flow models the LWR model is based on the mass conservation, 
i.e., traffic conservation, and is described by a first-order, nonlinear PDE: 

Pt + f(ph = 0, (i.i) 

in which f(p) = pt>*(p) is the traffic flow rate. Equation (1.1) is in conservation form. 
It is also called a kinematic wave model since it shows wave properties analogous to 
those of gases. There are many numerical methods to solve the LWR model. One 
approach is to solve the Riemann problem and apply a Godunov method for this 
model. Both solutions to the Riemann problem and the Godunov method are well- 
developed for hyperbolic conservation laws (Smoller, 1983). Another approach is to 
use the demand and supply functions (Lebacque, 1996; Daganzo, 1995), which turns 
out to be variants of Godunov's method. 

The PW model, derived based on microscopic car-following models, discards the 
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equilibrium assumption. It is a second-order system of hyperbolic conservation laws 
with a source term: 

Pt + (pv) x = 0, (1.2) 

v t + vv x + C -°p x = V APLZ^ ) (1.3) 

P T 

in which the constant cq is the traffic sound speed and r is the relaxation time. 
Equation (1.2) is the continuity equation, and (1.3) is the momentum equation. The 
PW model relates to driver behavior models better than the LWR model because 
it accounts for drivers' anticipation and inertia. However it's been shown that the 
LWR model is an asymptotic approximation of the PW model (Schochet, 1988). The 
PW model better captures non-equilibrium wave phenomena in traffic flow. Another 
property of the PW model is that it is unstable under certain situations. Since there 
are no known analytical solutions to the PW model, we use numerical methods to 
solve it in this research. All of these methods are developed from Godunov's method. 
Also based on microscopic models, Zhang (1998) developed a new non-equilibrium 
traffic flow theory : 

Pt + (pv) x = 0, (1.4) 

{pviip)) 2 v*(p) -v 

v t + vv x ^ p x = . (1.5) 

P r 

In the momentum equation (1.5), a varying sound speed c = pv^(p) has been in- 
troduced. It has been shown that this new model avoids "wrong-way travel" which 
is exhibited in the PW model (Zhang, 1998), and the model is always stable. Wave 
solutions to this model were discussed in (Zhang, 1999a), and finite difference approx- 
imations were studied in (Zhang, 2000a). In this research we perform the numerical 
simulations of this model. 

For a roadway with inhomogeneities such as a change in number of lanes, curva- 
ture and slopes, the LWR model can still be used, but the equilibrium speed-density 
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relationship varies with location. By introducing an inhomogeneity function a(x) 
which is a profile of the roadway at the location x, we can write the inhomogeneous 
LWR model as 

Pt + f(a,p) x = 0. (1.6) 

Here the traffic flow rate f(a,p) is a function of the inhomogeneity a(x). By writing 

a t = 0, (1.7) 

the inhomogeneous LWR model is a non-strictly hyperbolic system. There have 
been empirical methods to solve the inhomogeneous LWR model (Lebacque, 1995; 
Daganzo, 1995a). In this research, we develop a rigorous procedure to solve the LWR 
model based on the work by Isaacson et al. (1992) and Lin et al. (1995). We find 
the solutions that are consistent with those given by Lebacque. However, our method 
can be extended to solve higher-order inhomogeneous models while those of Lebacque 
and Daganzo cannot. 

Multi-commodity models are discussed in (Daganzo, 1994 and 1995; Jayakrish- 
nan, 1991; Vaughan et al., 1984). In these models traffic flow is disaggregated by 
origins, destinations or departure times. All of these models are based on traffic 
conservation. A First-In-First-Out (FIFO) discipline is assumed in all of these multi- 
commodity models. The model studied by Vaughan et al. is a continuous model. 
The models studied by Jayakrishnan and Daganzo are discrete models. In these two 
discrete models vehicles that are close to each other (in the sense of location or time) 
and have the same origin or destination or some other common characteristics are 
combined as a macroparticle. The macroparticles in a zone are ordered by time or 
location. Macroparticles are moved according to traffic conditions, which are solved 
with a link flow model. In the thesis, we introduce a new multi-commodity model 
that is more efficient in moving macroparticles. 
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1.2.2 Hyperbolic Systems of Conservation Laws and Godunov Methods 

A system of hyperbolic conservation laws (Smoller, 1983) takes the following form: 

u t + f(u) x = 0, (1.8) 

where u = (ui, ••-,«„) G R n ,n > 1, and (x, t) G R x R + . The n eigenvalues of 
the differential of f(u), Df(u), are denoted as Ai, • • • , A„. The solutions related to 
i-th eigenvalue are called i- family wave solutions. If the eigenvalues are distinct, the 
system (1.8) is a strictly hyperbolic system. To solve (1.8), initial and boundary 
conditions are needed. The Riemann problem is to solve (1.8) with the following 
jump initial condition: 

{Ui, x < 
(1.9) 
u r , x > 

where ui and u — r are constants. 

It is well-known that the weak solutions to the Riemann problem exist and are 
unique for the system (1.8) under the so-called "Lax's entropy condition" (Lax, 1972). 
The system admits discontinuous solutions, i.e., shock waves, and the wave speed s 
is determined by Rankine-Hugoniot condition: 

s[u] =[/(«)], (1.10) 

where [u] = ui—u r , and similarly, [/(«)] = f(ui)—f(u r ). For valid i-family shock wave 
solutions, the entropy inequalities hold, i.e., Aj(ti/) > s > \i(u r ). There are continuous 
solutions u = u(£,), £ = x/t, which satisfy the ordinary differential equations 

-£«* + /(«)* = o. (l.n) 

For most general initial and boundary conditions, there are no analytical solu- 
tions to (1.8), hence, one must use numerical methods to solve it. The Godunov 
method is one of the most efficient numerical methods, which combines solutions to 
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the Riemann problem and conservation laws. In a Godunov method, the space region 
[a, b] is divided into N grids x = a, xi, • • ■ ,xn-i,xn = b; the time scale [t ,ti] are 
partitioned into M time steps t° = to,t l , ■ ■ ■ ,t M ~ l ,t M = t±. The hyperbolic system 
of conservation laws (1.8) can be approximated by the finite difference equations: 

ur - u[ , m-1/2) - w+1,2) _ n ni2 . 

At + Ax ~ U ' [ } 

where U{ is the average of u(x, t) in grid [xi-1/2, Xi+1/2] at time step tj, similarly C// + 
is the average at time step tj+i, U*_ x , 2 is the average of u(x,t) during time interval 
[tj,tj + i] at grid boundary Xi-1/2, similarly U* +1 , 2 is the average at boundary x i+ i/ 2 - 
The boundary flux f{U*_ 1 i 2 ) is calculated by solving a Riemann problem at each cell 
edge with the following initial conditions: 

, Wv f U i-1, X < Xi-1/2 

[ U J i, X> X i+1/2 

Equations (1.12) says that the increment of u is equal to the difference between both 
boundary fluxes, which is the general idea of conservation. 

In this thesis our second-order models are not exact conservation laws since they 
have a source term, therefore our methods have been extended to address the effect 
of source terms. 

1.3 Layout of the Thesis 

In chapter 2, we study the homogeneous LWR model through theoretical discussions 
and numerical simulations. This is the first step for us to understand traffic flow 
models and associated numerical methods. In chapter 3, Godunov- type methods 
are developed for Zhang's model. The Riemann problem is solved in detail, and 
computational results are discussed. In chapter 4, the PW model is studied with 
several different Godunov-type finite difference methods, and the stability of the 
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model is tested. In chapter 5, the Riemann problem for the inhomogeneous LWR 
model is solved rigorously. In chapter 6, a multi-commodity model based on the LWR 
model is studied. In chapter 7, possible extensions of this research are discussed. 
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Chapter 2 

The LWR Model and Its 
Numerical Solutions 



The landmark paper by Lighthill and Whitham (1955) set the tone for many re- 
searchers' investigations into the theory of traffic flow. Introduced for traffic flows on 
a single, long, and rather idealized road, the LWR model proposes that their dynamics 
is described by the following PDE: 

Pt + f(p)x = 0, (2.1) 

where subscript t means the partial derivative with respect to time t, and subscript 
x means the partial derivative with respect to location x. In (2.1), the function 
f(p) = pv*(p) is called the fundamental diagram of traffic flow, in which v*(p) reflects 
the equilibrium relationship between travel speed and traffic density. It's generally 
assumed that f(p) is a concave function, i.e., f pp {p) < 0. The characteristic wave 
speed A(p) = f p (p) = f*(p) + pv'*(p)- A(p) can be positive or negative. 

Weak solutions (Smoller, 1983) for (2.1) satisfy the integral form of the conser- 
vation law: 

^ p p(x,t)dx = f(x u t)-f(x 2 ,t). (2.2) 
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First, we discuss the Riemann problem for (2.1) theoretically and then give the nu- 
merical solutions. 

2.1 The Riemann problem 

We consider the Riemann problem for the LWR model with the following jump initial 
condition: 

{07, x < 
(2.3) 
p r , x > 

There are two types of wave solutions to the Riemann problem. The discontinuous 
solution is a shock wave 

{pi, x/t < s 
(2.4) 
p r , x/t > S 

where s is the shock wave speed. The shock wave speed is determined by the Rankine- 
Hugoniot jump condition, 

s = im. „ 5) 

The wave speed of a valid shock wave solution has to satisfy the entropy condition: 

X( Pl ) > s > X(p r ). (2.6) 

Specifically, for a concave fundamental diagram, a shock wave is a solution to the Rie- 
mann problem for (2.1) with initial conditions (2.3) when p\ < p r ; i.e., the upstream 
traffic density is lower. 

When the upstream traffic density is higher, solution to the LWR model with 
initial data (2.3) is a continuous rarefaction wave. The rarefaction wave is given by 
P = p(0>£ = x /^i where p satisfies the ordinary differential equation 

-fo + /(*>)€ = 0. (2.7) 
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If p(£) 7^ 0, we obtain 

A(p(0) = e, (2.8) 

from which we can find p(£) uniquely. On any characteristic x/t = £, p is constant. 

2.2 Computation of boundary fluxes 

Given a wave solution to the Riemann problem, we can compute the average p* at 
x = and therefore the flux f(p*) through the boundary. There are the following 
five cases: 

Case 1 When the solution to the Riemann problem is a shock wave with wave speed 
s > 0, p* = p t . 

Case 2 When the solution to the Riemann problem is a shock wave with wave speed 

S < 0, p* = p r - 

Case 3 When the solution to the Riemann problem is a rarefaction wave, and X(pi) > 0, 

P* = Pi- 
Case 4 When the solution to the Riemann problem is a rarefaction wave, and X(p r ) < 0, 

P* = Pr- 

Case 5 When the solution to the Riemann problem is a rarefaction wave, X(pi) < and 
X(p r ) > 0, p* is the solution of X(p*) = 0. 

Given initial and boundary conditions, we can use a first-order Godunov method 
to calculate the traffic conditions. The numerical solutions are given in next section. 
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2.3 Numerical solutions to the LWR model 

In this section, we solve the Riemann problem for (1.8) numerically. We use Newell's 
equilibrium model, 

v*(p) = Vf(l-exp{^(l- Pj /p}j (2.9) 

Without loss of generality, we set Vf — 1, c 3 ■, — 1, p 3 ■, — 1 to obtain 

v„(p) = l-exp(l--) (2.10) 

P 

The valid range for p and v is < p, v < 1. The corresponding flow rate, /* = pt> *, is 
a normalized fundamental diagram. 

1. Given the initial conditions: 

f 0.65 xe [01,2001] 
p(x,0) = \ (2.11) 

I 0.4 otherwise 
v(x,0) = v,(p(x,0)), VxE [01,8001] (2.12) 

we obtain the solutions shown in Figure 2.1, Figure 2.2. 

2. Given the initial conditions: 

f 0.65 xe [01,2001] 
p{x,0) = I (2.13) 

J 0.9 otherwise 
v(x,0) = v*(p(x,0)), \/xe [01,8001] (2.14) 

we obtain the solutions shown in Figure 2. 3, Figure 2.4. 

From these numerical solutions we can see the formation of shock waves and 
rarefaction waves. These solutions are only an approximation of real solutions. For 
example, the solutions after t — in Figure 2.4 are not exact jumps, while the 
theoretical solution to the LWR model with initial conditions (2.3) is still a jump at 
any time. However as Ax — > and At — > 0, the solutions given by Godunov's method 
converge to the exact solutions. 
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R1 : evolution of density p for LWR model 
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Figure 2.1: Rarefaction wave solution of the LWR model with initial data (2.3) 



solutions of p at different times 




Figure 2.2: Solutions from Figure 2.1 at selected times 
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H1 : evolution of density p for LWR model 




H1 : evolution of velocity v for LWR model 




100 200 300 



400 500 600 700 

x/l 



0.4 

0.35 

0.3 

0.25 

0.2 

0.15 



Figure 2.3: Shock wave solution of the LWR model with initial data (2.3) 
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Figure 2.4: Solutions from Figure 2.3 at selected times 



Chapter 3 

Zhang's Second-Order Traffic Flow 
Model and Its Numerical Solutions 



3.1 Introduction 

A theory of non-equilibrium traffic flow has been developed by Zhang (1998), and 
in Zhang (2000a) he developed the Godunov-type finite difference equations (FDE) 
for this model. Zhang's model is a second-order model, and can be written in the 
conservation form 

where <f>(p) is a velocity flux function and defined as 

<P'(p) = — =P«(p)) 2 - (3.2) 

P 

Here c(p) = —pv'^(p) is the traffic sound speed. 

In (3.1), v*(p) is the equilibrium speed. Some often-used equilibrium travel speed 
functions are listed along with the corresponding velocity flux functions <f>(p) as fol- 
lows: 
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Functions 


v*(p) 


<t>{p) 


Greenshields 


M l ~ PlPi) 




Polynomial 


v f (l-(p/ Pj r),n>l 


V i{PlP 3 f n 


Greenberg 


v ]n(pj/p) 


v 2 ln(p) 


Underwood 


v f exp{-p/pj) 


-vj(l + p/pj)exp{-p/Pj) 


Newell 


v f [l-exp(^(l- Pj /p))] 


?(^J + |)exp(2M(i- p ./ p) ) 



The equilibrium travel speed i>* is decreasing with respect to traffic density; i.e., 
vj{p) < 0. The fundamental diagram /*(p) = pv*(p) is concave; i.e., f'J(p) < 0. In 
(3.1), r is the relaxation time and the relaxation term v * ( - p >~~ v constrains the differ- 
ence between the real travel speed v and the equilibrium travel speed v*. When the 
relaxation term is 0, Zhang's model reduces to the LWR model: 

Pt + (pv*)x = 0. (3.3) 

Zhang's model has three different wave velocities (relative to the road): A first- 
order wave velocity and two second-order wave velocities. The first-order wave veloc- 
ity is the wave velocity of the corresponding LWR model: 

K(p) = v*(p) - c(p) = v*(p) + pvj(p). (3.4) 

The two second-order wave velocities are 

Ai,2(p,u) = vtc(p) = v±pvj{p). (3.5) 

The relationship between these three wave speeds along (p, v«) phase curves is that 

Ai = A* < A 2 . (3.6) 

The waves with wave speed Ai are called 1-waves; similarly the waves with wave speed 
A2 is called 2-waves. Since p, v > and v'^ < 0, the 2-wave speed A2 > for any p, v. 
Since Zhang's model is a hyperbolic system of conservation laws with a relaxation 
term, the system is stable when (Liu, 1979, 1987; Chen et al., 1994) 

Ai < A* < A2- 
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This condition is satisfied by Zhang's model. Thus Zhang's model is always stable. 

In the following sections we study Godunov-type methods and use them to solve 
(3.1) numerically. In section 2 we discuss Godunov's method and properties of Zhang's 
model. In section 3 we present a second-order Godunov method. In Section 4 we solve 
the Riemann problems numerically and discuss the order of accuracy for different 
methods. 

3.2 Godunov's method 

A Godunov-type finite difference method for Zhang's model was first presented in 
(Zhang, 2000a). In this section we review this Godunov method and solve the asso- 
ciated Riemann problem. 

The Godunov-type FDEs for Zhang's model are 

7+1 7 *3 *3 *j *3 

Pi ~Pi Pi+l/2 V i+l/2 ~ Pj-l/2 V i~l/2 , v 

k + h " ' [ } 



vj - vj 2 + yyPi+i^) 2 vyPi-1/2) 

k h 



In these FDEs, p\ is the average of p in cell i at time step j; i.e., 



(3.J 



p\ = - p(x,tj)dx. (3.9) 

Similarly vj is the average of v. We use p* 3 +1 / 2 as the average of p through the cell 
boundary Xi+1/2 over the time interval (tj,tj+i), i.e., 



P*U/2 = T f p(x i+1/2 ,t)dt. (3.10) 



Similarly we define v*^_ 1/2 , p* J _ 1 , 2 ,v*^_ 1 , 2 as boundary averages. 
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By measuring the source term with values at time t,+i, we write the evolution 
equations for Zhang's model as 

k ■ ■ 

pi +1 = p\- i(p?+i/2 v *ii/2 - p*-i/2<-i /2 ) ( 3 - n ) 

/ * J \ 2 / *7 \ 2 

j+1 1 r j kA V i+l/2) ,, *j n ^i-l/2J ,, *j u 

v i = JY+^) { l ~ h^ 2 + ^^+1/2) 2 ^Pi-1/2)] 

+-v.(A +1 )} (3-12) 

T 

Provided traffic conditions (p, f) at time t,-, traffic conditions at time £,- + i can be 
calculated if we know the boundary averages P*+i/2i v *+i/2i Pi-i/2i v i-i/2- The com P u_ 
tation of P*+i/2i v *+i/2 a ^ the cen boundary £1+1/2 during the time interval (tj,tj+i) 
depends on a Riemann problem for (3.1) with the following initial conditions 

[/,, if x - x i+1/2 < 
iii+i/2(a:,tj) = ^ , (3.13) 

C/ r , if x - x i+1 / 2 > 

where we define the state variable u(x,t) = (p,v), U- = (pl,vl) = u(xi,tj) and left 
and right states U\ = (pi,vi), U r = (p r ,v r ). In a first-order Godunov method, we use 
the cell averages Pi — p\,vi — v\ as the left side (upstream) traffic conditions and 
p r = pi+i,v r = vf +1 as the right side (downstream) conditions. In a second-order 
Godunov method, we use higher-order approximations to the left and right states. 
(For details for a second-order Godunov method, refer to Section 3.3.) 

Here we have neglected the relaxation term in (3.1) when solving the Riemann 
problem. This Riemann problem has been discussed by Zhang (1999a), and the 
solutions to the boundary averages are provided there. The solutions are self-similar 
and can be expressed in the form of 
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3.2.1 Solutions of the boundary averages 

There are 8 types of wave solutions to the Riemann problem, which are combinations 
of two 1-waves and two 2-waves. The calculation of the boundary averages depend on 
the type of solutions. The formula for calculating the boundary averages are listed 
as follows. 



1. The wave solution is a 1-shock when the initial conditions satisfy 

HI: v r — Vi = ■ 
The wave speed is 



2 (P* - Pr)(<P(pl) -0(Pr)) 
Pi + Pr 



p r > pi, V r < V t . (3.14) 



p r v r - PiVi 

S = 

Pr ~ Pi 

The boundary averages {p*i + i/2i v l+i/2) are gi ven m the following table: 



(3.15) 



HI 


„ _ p r v r -pivi 

Pr-Pl 


Pi+l/2 


*j 

V i+l/2 


s> 


Pi 


Vl 


s< 


Pr 


V r 


s = 


Pl+Pr 

2 


Vl+V r 

2 



2. The wave solution is a 2-shock when the initial states satisfy 

H2: v r — Vi = 
The wave speed is 



/ 2 (p / -p r )(0(p / )-0( Pr )) 

Pi + Pr 



p r < pi, V r < Vi (3.16) 



p r v r - pivi 



>0. 



The boundary averages (p*+ 1 / 2 5 ,y i4.i/2 



Pr - Pi 

) are given in the following table: 



(3.17) 



H2 


„ _ PrVr-pl'Oi 
Pr-pi 


*j 

Pi+l/2 


V i+l/2 


s> 


Pi 


Vl 
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3. The wave solution is a 1-rarefaction when the initial states satisfy 



Rl: v r -vi = v*(p r ) - v*(pi), p r < Pi, v r > Vi 



The characteristic speed of a 1-rarefaction wave is 



Xi(p,v) = v + pvj(p) 



(3.18) 



(3.19) 



The boundary averages are the left state when Ai(p;,i>;) > 0, similarly they are 
the right state when Ai(p r ,t> r ) < 0. Otherwise, {p* 3 +1 /2i v *i+i/2) are ^ e solutions 
of the equations 



A i(Pm/ 2 . <+i/ 2 ) = P*Ii/2 v *'(pZi/2) + <+i/ 2 = ° 

v *ii/2- v i = v Ml 1/2 )- v *(pi)- 

We simplify equations (3.20,3.21) as 

A *(P*+i/ 2 ) = v *{Pi) ~vi = Av 

1/2 = MP*U/ 2 )- Av - 



J i+1/ 



The boundary averages (p*I 1 / 2 ) t 'i4-i/2) are gi ven m the following table: 



(3.20) 
(3.21) 

(3.22) 
(3.23) 



1+1/2' "1+1/2 



Rl 


Ai 


Pi+l/2 


V i+l/2 


^i(pi,vi) > 


Pi 


Vl 


Xl(p r ,V r ) < 


Pr 


v r 


o.w. 


solution to (3.22,3.23) 



4. The wave solution is a 2-rarefaction when the initial states satisfy 



R2: v r - vi = v*(p t ) - v*(p r ), p r > p h v r > v t 



The characteristic speed of the 2-rarefaction wave is 



(3.24) 



M(p,v) = v - pvj(p) > 0. 



(3.25) 
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The solutions of (/0* + i/2j u *+i/2) are gi ven i n the following table: 



R2 


A 2 


*j 
Pi+l/2 


*j 

V i+l/2 


A 2 >0 


Pi 


Vl 



5. The wave solution is a 1-rarefaction + 2-rarefaction when there exists an inter- 
mediate state (p m ,v m ) satisfying 

Rl: v m - vi = v*(p m ) - v*(pi), p m < pi, v m > vi (3.26) 

R2: v r -v m = v*(p m ) - v*(p r ), p r > p m , v r > v m . (3.27) 

That is to say, p m satisfies 



2*f*(p m ) -V*(pi) - V*(p r ) - (v r -Vi) = 



(3.28) 



in which p m < pi, p m < p r . We can write v m as 



V m =V*(p m ) +Vi -V^(pi). 



(3.29) 



The boundary averages {p* i 3 + i/ 2 -i v *i+ii2) are given in the following table: 



1+1/2' u i+l/2 



R1-R2 


Ai 


Pi+l/2 


*j 

V i+l/2 


^i(pi,vi) > 


Pi 


Vl 


Al(p m ,W m ) < 


Pm 


Vm 


o.w. 


solution to (3.22,3.23) 



6. The wave solution is a 1-rarefaction + 2-shock when there exists an intermediate 
state (p m , v m ) satisfying 

Rl: v m - vi = v*(p m ) - v*(pi) ,p m < pi, v m > v { (3.30) 



m-.v r -v m = - ^(^-Prm^^ Ml , Pr<Pm , Vr<Vm . (3.31) 
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That is to say, p m satisfies 



2(pm-Pr)(4>(pm)-4>(Pr)) 
Pm+pr 



- ( Vr - Vl ) =0 (3.32) 



in which p r < p m < p\. We can write v m as 

V m =V*(p m ) +Vl ~V*(pi). 

The boundary averages (p* J +1 /2, 1 C+1/2) are gi ven as the following 



(3.33) 



R1-H2 


Ai 


Pi+l/2 


V i+l/2 


Xi(pi,vi) > 


Pi 


Vl 


><l(Pm,V m ) < 


Pm 


Vm 


o.w. 


solution to (3.22,3.23) 



7. The wave solution is a 1-shock + 2-shock when there exists an intermediate state 
(Pm,v m ) satisfying 

m:v m -v t = - ^-^)«M- &S ,p m> p hVm <v l (3.34) 
m:v r -v m = - J ^-^)S ,A-<P m ,«r<« m - (3-35) 



That is to say, p m satisfies 



2{pl-p m ){4>(pi)-4>{pm)) __ I 2{pm-pr){<t>{pm)-<j>{Pr)) 



Pl+pm V pm+pr 

in which p m > Pi, Pm > Pr- We can compute t> m as 



- (^r - ^ 



2(pz - Pm)(4>(pi) ~ 4>(Pm)) 



Pi + Pr, 



+ ^- 



0(3.36) 



(3.37) 



The boundary averages (p*+i/ 2 ,v*+i/ 2 ) are given in the following table: 



H1-H2 


„ _ PmV m -plVl 
Pr-pl 




*j 

V i+l/2 


s> 


Pi 


Vl 


s< 


Pm 


Vm 


s = 


Pl+Pm 
2 


Vl+V m 

2 
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The wave solution is a 1-shock + 2-rarefaction when there exists an intermediate 
state (p m , v m ) satisfying 



m-.v m -v l = - ^fo-^Wta)- feir ,p m >p l ,v m <v l (3.38) 

V*{Pm) ~ V*{Pr) ,Pr > Pm, V r > V m . (3.39) 



R2: v r — v m 



That is to say, p m satisfies 

in which p m > p h p m > p r . We compute v m as 



(3.40) 



2(pi - p m )(4>(pi) -<f)(p m )) 



Pi + Pm 



+ Vi. 



(3.41) 



The boundary averages (p* J +1 / 2 , v iii /o) are gi ven m the following table: 



'i+l/2' u i+l/2 



H1-R2 


„ _ PmVm-plVl 
Pr-pl 


*j 

Pi+l/2 


V i+l/2 


S > 


Pi 


Vl 


s < 


Pm 


Vm 


s = 


pl+Pm 
2 


Vl+Vm 
2 



3.2.2 Some points concerning the implementation of the Godunov method 

In subsection 3.2.1 above, we studied the solutions of the boundary averages. We 
now discuss how to get a stable, convergent and efficient numerical method. 

In a linear hyperbolic system, Godunov's method is stable and convergent if the 
Courant-Friedrichs-Lewy (CFL) number is less than unity. Similarly we require the 
CFL number be less than unity for Zhang's model; i.e., 

k 



max 



h 



X 2 (p,v) 



< 1 



(3.42) 



where \2{p, v) — v — pvj(p) has the bigger magnitude of two wave velocities. 
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There are two types of boundary conditions (BC) : Dirichlet BC and Neumann 
(or natural) BC. If the interval for x is [a,b], we will impose boundary condition on 
u(a — |, tj) and u(b + |) instead of the real boundary x = a and x = b. This is to say 
we have to solve a Riemann problem to get the fluxes on both the end boundaries 
instead of imposing boundary conditions directly on those fluxes. 1 

The source term, s(u) = (0, v * l - p >~ v ^ doesn't involve spatial gradients and there- 
fore remains bounded as terms on the left hand-side of (3.1) go to oo. We approximate 

■_l_1 j j 3 ' + 1 _|_ j j 3 

the source term implicitly with s(U? ) or s(— *— ^ — L ) in order to improve the stability 
property of the Godunov method. However, when the source term is stiff; i.e., when 
r is small, the problem of numerical instability may arise. In this case we use a much 
smaller time increment k <^ t. 

The solution to the Riemann problem is important both theoretically and com- 
putationally. One can compute the numerical solutions when all of the Riemann 
problems are well-posed and solvable at each step of the iteration. However, when 
the left and right states for a Riemann problem are far from each other, the inter- 
mediate state (p m ,v m ) for the wave solutions may be out of domain of validity, e.g., 
p m < 0. In this case, we have a "vacuum problem" and hence the numerical solutions 
can not be uniquely determined. 

The cost of solving the Riemann problem determines the computational efficiency 
of the numerical method. Here we propose improvements to the computational ef- 
ficiency for Zhang's model. In Godunov's methods for Zhang's model, most of the 
calculations are basic arithmetic computations except the calculation of the interme- 
diate state (p m ,v m ) in equations (3. 28), (3. 32), (3. 36) and (3.40) and the solutions of 
equations (3.22) and (3.23)). Given (pi,vi) and (p r ,v r ), the nonlinear algebraic equa- 
tions (3. 28), (3. 32), (3. 36) and (3.40) can all be written in the form of g(p m ) = 0. The 
functions g{p m ) for these equations are all monotonically decreasing in the interval of 

1 For the detailed discussions on treatment of boundary conditions, refer to (Zhang, 2000a). 
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validity of p m . To find p* 3 +1 / 2 from equation (3.22), we define a function 

g (p) = A*(p) - Av, p r {p m for R1-R2,R1-H2) < p < Pl . (3.43) 

Since /*(p) is concave and A*(p) = fl(p), we find A*(p) and g(p) are decreasing. Here 
we chose secant method to solve these equations since it is very efficient when the 
associated functions are monotonically decreasing. 

3.3 A Second-order Godunov Method 

In this section we introduce a second-order finite difference method for Zhang's model. 
This method is a two-stage predictor/corrector method. In this method, (3.1) is 
decoupled into two nonlinear scalar equations in each interval (2^-1/2,3^+1/2) & t time 
tj and the predictor/corrector procedures are applied to those scalar functions. 
We can write Zhang's model as: 

u t + A{u)u x = s(u) (3.44) 

where 

I P( X >*) I /O ACS 

u = (3.45) 

\ v(x,t) J 

and 

A(u) = ( V P I (3.46) 

y pvj 2 v J 

The two eigenvalues and corresponding eigenvectors are 

Ai(w) = v + pv*(p), r x (u) = [l,vj(p)]\ 
\ 2 (u) =v- pvj(p), r 2 (u) = [1, -vj(p)} 1 . 
We diagonalize A{u) by 



(3.47) 



. Ai(m) 
T-\u)A{u)T{u) = I V ' |=A(w), (3.48) 

A 2 (w) 
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where the transformation matrix T(u) is 

T{u) = ( l I. (3.49) 

\ vj(p) -vj(p) J 

Letting W = T _1 (w)-u, Zhang's model (3.1) under the transformation becomes 

W t + K{u)W x = r _1 («)5(«). (3.50) 

For the solution w(x,t) to a scalar equation Wt + X(w)w x = 0, the first-order 
Godunov method uses a step function Wi(x,tj) to interpolate the solution, i.e., 

w I (x,t j ) = wf, if Xj_i/ 2 < x < Zj+i/2. (3.51) 

In a first-order Godunov method, the Riemann problem has the following jump initial 
conditions: 

W ^ /2 = Wl (3.52) 

w i-i/2 = <■ 

For a second-order Godunov method, we interpolate the initial condition with a piece- 
wise linear function 

Wi(x,tj) = w{ -\ A VL wf, if Xi-1/2 < x < Xi+i/2- (3.53) 

Then we do a half-step prediction 

<^' L = ^ + |(i-ak)!)a^V (354) 

<-^ R = ^'-|(H-AK)|)A^X, 
where A yL u^ is the van Leer slope defined as (all subscripts j have been suppressed) 

A VL w, = 



Si ■mm(2\w i+1 - Wi\, 2\w t - Wj_i|, ^\w i+1 - Wi-i\), £ > 



0, otherwise 

Si = sign(w i+ i - Wi-i) (3.55) 

f = (Wi+i - Wi) ' (wt - Wi-i) (3.56) 



CHAPTER 3. Zhang's Second-Order Traffic Flow Model and Its Numerical Solutions 27 

The van Leer slope limiter ensures that the method remains second order when the 
solution w(x,t) is smooth and eliminates Gibb's phenomenon at discontinuities. 

We apply the procedure above to the two scalar equations of the related homo- 
geneous 2x2 system in (3.50) to obtain W? +1 L ' and W^L ' . Given the half-step 
values of W- +1 L ' and W/_ 1 L ' , U 3 i+l L ' and U^L ' can be calculated by an inverse 
transformation: 



u^ L = T{urm:T 

T j+1/2,R _ , j, j+l/2,R 



(3.57) 



We then solve the Riemann problem to find the boundary averages. 

3.4 Numerical Solutions of Zhang's model 

Based on the discussions in the former sections, we carry out some numerical com- 
putations to test the validity of the Godunov method and the properties of Zhang's 
model. 

Here we use Newell's equilibrium model, 

<p) = t;/(l-exp{M(l-p,/p)}). (3.58) 

and set Vf = l,Cj = l,pj = 1 to get the standardized equilibrium relationship: 

v*(p) = l-exp(l--) (3.59) 

P 

The domain for pis0<p< l 2 , the range for v* is the same. The traffic flow rate is 
/* = pv*. The equilibrium travel speed i>*(p) and flow rate /*(p) are shown in Figure 
3.1. 

The subcharacteristic; i.e., the first-order characteristic velocity, (i.e., the wave 
velocity of the corresponding LWR model) is 

A, = 1- (l + -)exp(l- -), 
P P 



J We use the limit of i>* when p — > as its value at p = 0. 
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and the eigenvalues are 

Ai )2 = u±-exp(l ). 

P P 

When < p < 1, we have 

|A 2 | < \v\ + 1. 

The CFL condition number is defined as: 

max |— A 2 (p, v )| < — (maxt> + 1). (3.60) 

Since the CFL number is no larger than 1, we find 

k < . (3.61) 

max v + 1 

In all of the computations that follow, we let x G [0/, 800/], where / is the unit of 

length. Here the number of grid points is denoted by N, and h = ^^ is the space 

step. We let T = Kr denote the final time (r is the unit of time), and m the number 

of time steps and k — — . From the CFL condition, we have 

(maxv + l)KNr 

± ^ < 1. 3.62 

800 m I ~ 

Setting t = I = 10.0, and m — N, CFL condition is not violated when K < 400 since 
maxi) < 1. 

3.4.1 Riemann solutions 

In the following four computations we use the first-order Godunov method to examine 
different types of waves solutions. With four well-chosen jump initial conditions, we 
can observe four different type of waves H1-H2, R1-R2, R1-H2 and H1-R2. To prevent 
the 2-waves from relaxing to 1-waves in a short time, here we rescale the relaxation 
time t — ► lOOOr. For each computation we present a contour plot of both p and v 
until To = 400r, and several 2-D curves at different time t have been selected from 
the contour plot. 
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density - velocity relationship 
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Figure 3.1: Newell 's Equilibrium p-v^j p-f* Relationship 
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Computation 1 We use the following initial conditions: 

p(x,0) = 0.65, VxG[0Z,800Z] (3.63) 

f u,(0.65) xe [01,2001] 

v(x, 0) = < (3.64) 

I i>*(0.65) — 0.2 otherwise 

The solutions are of H1-H2 type, shown in Figures 3.2 and 3.3. These figures 
show that the downstream travel speed keeps increasing along with the prop- 
agation of the 2-shock wave. This is due to the effect of the relaxation term. 
We can predict that when the downstream travel speed reaches i>*(0.65), which 
is the equilibrium travel speed with respect to the downstream traffic density 
p = 0.65, the 2-shock disappears and a 1-rarefaction wave forms. 

Computation 2 We use the following initial conditions: 

p(x,0) = 0.65, VxE [01,8001] (3.65) 

f t,,(0.65) xe [01,2001] 

v(x,0) = < (3.66) 

I f* (0.65) + 0.2 otherwise 

The solutions are of R1-R2 type, shown in Figures 3.4 and 3.5. These figures 
show that the downstream travel speed keeps decreasing along with the propa- 
gation of the 2-rarefaction wave. This is also due to the effect of the relaxation 
term, which will relax the 2-rarefaction wave to a 1-shock wave. 

Computation 3 We use the following initial conditions: 

f 0.65 x E [01, 2001] 
p(x,0) = \ (3.67) 

I 0.4 otherwise 
v(x,0) = u„(0.65), \/x E [01,8001] (3.68) 

The solutions are of R1-H2 type, shown in Figures 3.6 and 3.7. Here the effect 
of the relaxation term will relax the 2-shock wave to a 1-rarefaction wave. 
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Computation 4 We use the following initial conditions: 

, 0.65 x E 101,2001} 
p{x,0) = { (3.69) 

0.9 otherwise 
v{x,0) = u,(0.65), Vx e [01, 8001] (3.70) 

The the solutions are of H1-R2 type, shown in Figures 3.8 and 3.9. Here the 
effect of the relaxation term will relax the 2-rarefaction wave to a 1-shock wave. 



3.4.2 A General Solution and Convergence Rates 

In this subsection the relaxation time is no longer rescaled since we wish to observe 
the solutions to (3.1) for a longer time. Here we are interested in the solution for 
< t < T = 400r. 

Using the general initial condition 

p(x,0) = 0.65 + sin (—)/4, (3.71) 

v(x,0) = v*(p(x,0)) + 0.1, \/xe [01,8001], (3.72) 

we use the first-order Godunov method to obtain the solutions shown in Figures 3.10 
and 3.11. From these figures we see a shock wave forms at the downstream section, 
and a complicated combination of rarefaction waves forms in the upstream section. 
The solutions also show that p-v are in equilibrium by t — lOOr, i.e., v = f*(p), due to 
the effect of the relaxation term, although the initial condition is not in equilibrium. 
Next, we compute the convergence rate for the first- and second-order methods 
with Neumann boundary conditions with initial conditions (3.71,3.72). The conver- 
gence rate is calculated from the comparison between the solution on different grids. 
The grid numbers 2N and N generate different grid sizes: | and h. We denote the 
solutions at time T as (L^ 2Ar )^ and (U^)^ respectively, and define the difference 
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H1-H2: evolution of density p 
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Figure 3.2: Left-Shock-Right-Shock Wave Solution 
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Figure 3.3: H1-H2 Solutions from at Selected Times 
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R1 -R2: evolution of density p 
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Figure 3.4: Lcft-Rarefaction-Right-Rarefaction Wave Solution 
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Figure 3.5: R1-R2 Solutions at Selected Times 
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R1 -H2: evolution of density p 
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Figure 3.6: Left-Rarefaction-Right-Shock Wave Solution 
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Figure 3.7: R1-H2 Solutions at Selected Times 
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H1-R2: evolution of density p 




Figure 3.8: Left-Shock-Right-Rarefaction Wave Solution 
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Figure 3.9: H1-R2 Solutions at Selected Times 
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vector (e 2Ar_7V )^ 1 between these two solutions as 

e™~ N = \(U™ 1 + Ul N )-U t N ^ = l,...,N. (3.73) 

Then the relative error between the two solutions is defined as the norm of the dif- 
ference vector: 

e 2N-N = || e 2M ||. (3.74) 

The convergence rate is defined as 

e 2N-N 

r = lo S2(^z^)- (3-75) 

In (3.74), the norm can be L 1 -, L 2 - or L°°-norm. 

For N equal to 64, 128, 256, 512 and 1024, the relative errors and convergence 
rates for the first-order method are given in Table 3.1, and those for the second- 
order Godunov's method are computed and given in Table 3.2. For the first-order 
Godunov's method, the convergence rates related to L 1 -norm errors are around 1, 
which is larger than that related to L 2 -norm errors and even larger than that related 
L°°-norm errors. From Table 3.1 we also see that the rates for p and v are consistent 
since p and v are in equilibrium at 400r. In Table 3.2, the convergence rates for the 
second-order method, although better slightly than those for first-order method, are 
not second order due to the effect of the relaxation term. 
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A general solution: evolution of density p 




Figure 3.10: A General Solution of Zhang's Model with Neumann BC 
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Figure 3.11: The Solutions from Figure 3.10 at Selected Times 
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p 


128-64 


Rate 


256-128 


Rate 


512-256 


Rate 


1024-512 


L 1 


3.43e-03 


0.969 


1.75e-03 


0.985 


8.84e-04 


0.994 


4.44e-04 


L 2 


6.94e-03 


0.632 


4.48e-03 


0.635 


2.88e-03 


0.718 


1.75e-03 


L°° 


4.56e-02 


0.0219 


4.49e-02 


0.0973 


4.20e-02 


0.355 


3.28e-02 


V 


128-64 


Rate 


256-128 


Rate 


512-256 


Rate 


1024-512 


L 1 


4.83e-03 


0.973 


2.46e-03 


0.985 


1.24e-03 


0.991 


6.25e-04 


L 2 


9.86e-03 


0.631 


6.37e-03 


0.628 


4.12e-03 


0.695 


2.54e-03 


L°° 


6.53e-02 


0.0243 


6.42e-02 


0.0903 


6.03e-02 


0.313 


4.86e-02 



Tabic 3.1: Convergence rates for the first-order method with initial conditions (3.71,3.72) 



p 


128-64 


Rate 


256-128 


Rate 


512-256 


Rate 


1024-512 


L 1 


5.81e-03 


0.966 


2.98e-03 


1.03 


1.46e-03 


1.05 


7.06e-04 


L 2 


9.85e-03 


0.823 


5.57e-03 


0.836 


3.12e-03 


0.885 


1.69e-03 


L°° 


3.73e-02 


0.215 


3.22e-02 


0.232 


2.74e-02 


0.616 


1.79e-02 


V 


128-64 


Rate 


256-128 


Rate 


512-256 


Rate 


1024-512 


L 1 


8.24e-03 


0.973 


4.20e-03 


1.03 


2.05e-03 


1.04 


9.93e-04 


L 2 


1.40e-02 


0.842 


7.80e-03 


0.865 


4.28e-03 


0.833 


2.40e-03 


L°° 


5.40e-02 


0.250 


4.54e-02 


0.354 


3.56e-02 


0.527 


2.47e-02 



Table 3.2: Convergence rates for the second-order method with initial conditions (3.71)ini.2 



Chapter 4 

The PW Model and Its Numerical 
Solutions 



4.1 Introduction 

The Payne- Whitham (PW) model, suggested by Payne (1971) and Whitham (1974), 
is one of the first non-equilibrium traffic models. It can be written as 

Pt + pv x + vp x = (4.1) 

v t + w x + C -°p x = OMzl, (4.2) 

P r 

where the traffic sound speed c > is constant, and i>*(p) is the relationship between 
velocity v and density p for equilibrium states. The fundamental diagram f*(p) = 
pv*(p), reflecting basic features of a roadway, is assumed to be known. Setting m = pv, 
(4.1) and (4.2) can be written as 

m J, + U+*J. = l £ w' <4 ' 3) 

which is in conservative form 

u t + f(u) x = s(u), (4.4) 

39 
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with state u = (p, m) of traffic flow. 

For Co = pvl(p), system (4.3) becomes Zhang's model, which is discussed in Chap- 
ter 3. Both PW model and Zhang's model are "systems of hyperbolic conservation 
laws with relaxation" in the sense of Whitham (1959; 1974) and Liu (1987). 

Schochet (1988) has shown that, as r — > 0, the system (4.3) admits the limit 



dx" 1 



Pt + (pv*) x = v^-^. (4.5) 



Note that (4.5) is the LWR model with viscous right handside. 

Generally, the equilibrium velocity is assumed to be decreasing with respect to 
density; i.e., vj{p) < 0; the equilibrium flow rate is concave; i.e., f'J(p) < 0. Equation 
(4.3) has three different wave velocities (relative to the road): A*, Ai and A2. The 
first-order wave speed A* is 

A*(p) = v*(p) +pvj(p). (4.6) 

A* can be positive or negative. Since the wave speed of the degenerate system (i.e., 
the LWR model), it is called a sub-characteristic. The second-order wave speeds, or 
"frozen characteristic speeds" (Pember, 1993), are 

Xi(p,v) = v-c , (4.7) 

\ 2 {p,v) = v + c . (4.8) 

Since v, cq > 0, Ai < A2 and v + cq > 0. 

Whitham (1959; 1974) showed that the stability condition for the linearized sys- 
tem with a relaxation term is 

Ai < A* < A 2 when t = 0. (4.9) 

Liu (1987) showed that if condition (4.9) is always satisfied, then the corresponding 
LWR model is stable under small perturbations and the time-asymptotic solutions of 
the system (4.3) are completed determined by the equilibrium LWR model. Chen, 
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Levermore, and Liu (1994) showed that if condition (4.9) is satisfied, then solutions 
of the system tend to solutions of the equilibrium equations as the relaxation time 
tends to zero. Besides (4.9), that Ai 7^ for all x and t can serve as another stability 
condition, since, otherwise, the standing wave f(u) x = s(u) may be singular and can't 
always be solved. 

This chapter is organized as follows. After discussing the boundary averages in 
section 2, we study several different Godunov methods for the PW model in section 
3. In section 4, we present numerical solutions to the PW model. 

4.2 Computation of the boundary averages of p and v 

To develop Godunov's methods for the PW model, we first partition a piece of road- 
way, e.g., an interval of [a, b], into N zones, and then approximate current traffic 
conditions p and v with certain types of functions. At each zone boundary, the av- 
erages of p and v over time have to be computed at every time step for computing p 
and v at next time step. 

For a system of conservation laws, we can compute the boundary averages by 
solving a Riemann problem, which has been discussed by Smoller (1983). For the 
homogeneous version of the PW model, we can adopt the solutions by Zhang (2000a). 
Zhang developed solutions to the Riemann problem for the homogeneous version of 
Zhang's model, which is similar to the PW model, and obtained 8 types of wave 
solutions and the formula for the boundary averages of p and v related to each type 
of solutions. However, the PW model has one relaxation term. The Riemann problem 
for (4.3) with a source term is still open. In this section, we still calculate the boundary 
averages based on the solutions to a Riemann problem for the homogeneous version 
of the PW model. 

Liu (1979) discussed the Cauchy problem for a hyperbolic system of conservation 
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laws with source terms. Inspired by his study, we present another approach of com- 
puting the boundary averages of p and v by solving a Cauchy problem for (4.3). This 
method is presented in the second part of this section. 

4.2.1 Computing the boundary averages from the Riemann problem 

In this subsection, we study the Riemann problem for the homogeneous version of 
(4.3) with the following jump initial conditions: 

{Ui, if x < Xi+i/o 
, (4.10) 

U r , if x > Xi+i/2 

where the left state U\ = (pi,vi) and the right state U r = (p r ,v r ) are constant. For 
computational purpose, we are interested in the averages of p and v at the boundary 
x = £j + i/2 over a time interval At, which are denoted by p*+u 2 an d v* 3 +1 , 2 ; i.e., we 
want to find 

At I /-At 



Pili/2 = & J P(* = 0>*)<lt and v*l 1/2 = —j v(x = 0,t)dt. (4.11) 

For the homogeneous PW model, the equations of the characteristic curves are 
written in terms of (p, v) instead of (p, to), and the left and right initial values for v are 
given as v i = — and v r = '—. Thus the boundary average of to, w*i 1/2 = P*ii/2 t 'j*+i/2' 
can be easily obtained once we computed P^/2 an d v * 3 + i/2- 

Determined by the relationship between left state (pi,vi) and right state (p r ,v r ), 
there are 8 types of wave solutions to the Riemann problem, including 4 first-order 
waves and 4 second-order waves. The four first-order wave solutions are a 1-shock, a 2- 
shock, a 1-rarefaction and a 2-rarefaction. The four second-order wave solutions are a 
H1-H2 (Left-Shock- Right-Shock) wave, a R1-R2 (Left-Rarefaction- Right-Rarefaction) 
wave, a R1-H2 (Left-Rarefaction-Right-Shock) wave and a H1-R2 (Left-Shock-Right- 
Rarefaction) wave. 
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In the remaining part of this subsection, we define the velocity flux function (j>{p) 
as <j)(p) = CqP, discuss the eight types of wave solutions to the Riemann problem one 
by one, and present a table containing the boundary averages p*+ 1/2 and v* 3 +1/2 for 
each case. 



1. The wave solution is a 1-shock if the left and right states satisfy 

HI: v r — Vi = 
The wave speed is 



2 (_Pj - Pr)(<P(pl) ~0(Pr)) 

PiPt 



p T > pi, v r < vi. (4.12) 



PrV r - PlVl 



(4.13) 



Pr ~ Pi 

The solutions (p*+i/2i v *+i/2) f° r case 1 are summarized in the following table: 



HI 


„ _ prVr-piVl 
Pr-Pl 


*j 

Pi+l/2 


*j 

V i+l/2 


s>0 


Pi 


Vl 


s < 


Pr 


V r 


s = 


Pl+Pr 
2 


Vl+V r 

2 



2. The wave solution is a 2-shock if the left and right states satisfy 

H2: v r — Vi = 
The wave speed is 



/ 2 (p / - Pr )(0( A )-0( Pr )) 

PlPr 



p r < Pi, V r < Vi (4.14) 



PrVr - PlVl 
Pr ~ Pi 



>0 



(4.15) 



The solutions (p. 



*3 *3 

i+1/2' u i+l/2 



for case 2 are summarized in the following table: 



H2 


„ _ p r v r -p l v l 

Pr-Pl 


Pi+l/2 


V i+l/2 


s> 


Pi 


Vl 
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3. The wave solution is a 1-rarefaction if the left and right states satisfy 
Rl: v r -vi = v*(p r ) - v*(pi) p r < pi, v r > vi 



(4.16) 



The characteristic velocity is determined by the first eigenvalue of the system: 



\i(p,v) 



V-Cq. 



(4.17) 



If Xi(pi,vi) > 0, the boundary averages p*+ 1/2 an d v *A.xj 2 are the left initial 
values for p and v , similarly, if Ai(p r ,t> r ) < 0, they are the right initial values. 
Otherwise, {p* 3 +1 /2-> v *i+i/2) are solutions to the equations: 



M{p2,^,v*A 



<+l/2 - CO = 0, 



7 i+l/2' ^+1/2/ 

v iii/2~ v i = v *(p*ii /2 )- v *(pi), 

which can be simplified as follows, 



V i+l/2 



Co, 



V *(p?+l/2) = C - V l + v *(pl)- 



(4.18) 
(4.19) 

(4.20) 
(4.21) 



The solutions (p* 3 +l 
lowing table: 



*] „*] 



/2> V i+l/2< 



to (4.20,4.21) for case 3 are summarized in the fol- 



Rl 


Ai 


Pi+1/2 


*j 

V i+l/2 


Xi(pi,vi) > 


Pi 


Vl 


Xl(Pr,V r ) < 


Pr 


v r 


O.W. 


solution to equations 4.20, 4.21 



4. The wave solution is a 2-rarefaction if the left and right states satisfy 
R2: v r -v t = v*(pi) - v*(p r ) p r > pi, v r > v t 
The characteristic velocity is the second eigenvalue of the system: 



(4.22) 



A 2 (p,-y) 



v + c > 0. 



(4.23) 
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The solutions (pl+i/ 2 , v i+i/2) f° r case 4 are summarized in the following table: 



i+1/2' "-1+1/2 



R2 


A 2 


*j 
Pi+l/2 


*j 

V i+l/2 


A 2 >0 


Pi 


Vl 



5. The wave solution is a 1-rarefaction + 2-rarefaction with an intermediate state 
(Pm, v m) if the left, right and intermediate states satisfy 



Rl: v m - vi = v*(p m ) - v*(pi) p m < Pi, v m > vf, 
R2: v r - v m = v*(p m ) - v^(p r ) p r > p m , v r > v m . 



Adding (4.24) to (4.25) we find 



2 * v*{p m ) - v^(pi) - v*{p r ) - (v r - vi) 



for p m < pi, p m < p r - Thus 



V m =V*(Pm) +Vi -V^(pi). 



(4.24) 
(4.25) 



(4.26) 



(4.27) 



The solutions (p* i+ i/2,v*l l / 2 ) f° r case 5 are summarized in the following table: 



R1-R2 


Ai 


Pi+l/2 


*j 

V i+l/2 


\i(pi,vi) > 


Pi 


Vl 


Xl(Pm,V m ) < 


Pm 


Vm 


o.w. 


solution to equations 4.20, 4.21 



6. The wave solution is a 1-rarefaction + 2-shock with an intermediate state (p m , v ri 
if the left, right and intermediate states satisfy 



Rl: v m -v h = v*(p m ) - v*(pi), p m < pi, v rn > vi (4.28) 

H2: v r — v m = - 



2(pm-Pr)(4>(pm)-4>(Pr)) 



PmPr 



, p r < p m , V r < V m . (4.29) 
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These two equations yield 



v*( Pm ) ~ v*(Pi) ~ ^ 2 (^-^)W^)-S i - ( Vr - Vl ) = 0, (4.30) 



for p r < p m < pi- Thus 



V m = V*(p m ) +Vl- V*(pi). 



(4.31) 



The solutions (p/Lrt^f+irt) f° r case 6 are summarized in the following table: 



y i+l/2> i+1/2 



R1-H2 


Ai 




V i+l/2 


Ai(pj,vj) > 


Pi 


Vl 


Ai(/o m ,u m ) < 


Pm 


Vm 


o.w. 


solution to equations 4.20, 4.21 



7. The wave solution is a 1-shock + 2-shock with an intermediate state (p m ,v m ) if 
the left, right and intermediate states satisfy 

m:v m - Vl = - /'""^y- ^, Pm>Pi,v m < Vl ; (4.32) 



H2: f r — v m 



2(pm-pr)(4>(pm)-4'(Pr)) 



PmPr 



, p r < p rn , v r < v m . (4.33) 



These two equations imply 



2(pi - Pm)((j)(pi) -4>{Pm)) . /2(p m -Pr)(0(Pm) -0(Pr)) 



PlPr, 



PmPr 



+ (v r -v l ) = 0,(4.34) 



for p m > pi,p m > Pr- Thus 



2(pi - Pm){4>{Pl) - 4>{Pm)) 



PlPr. 



Vl- 



(4.35) 



The solutions (p*+i/2>%*+i/2) f° r case ^ are summarized in the following table: 
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H1-H2 


„ _ PmVm-plVl 
Pr-pl 




V i+l/2 


s> 


Pi 


^i 


s< 


Pm 


<^m 


s = 


Pl+Pm 
2 


2 



The wave solution is a 1-shock + 2-rarefaction with an intermediate state (p m , v r , 
if the left, right and intermediate states satisfy 



HI: v m -vi 



R2: v r — v r 



2(pi-p m )(<l>(pi)-(j>(pm)) 



PlPm 



, Pm > Pi, v m < Vl] 



(4.36) 



V*{p m ) - V*(p r ), p r > p m , V r > V m . (4.37) 



These two equations imply 



_^ 2( Pi -p m )(0to)- ^JT + v ^ pm) _ Vm(j)r) _ {Vr _ Vl) - (4. 38) 



for p m > PhPm> Pr- Thus 



2 (p; - p m )((f>(pi) - 4>{Pm)) 

PlPm 



+ v L . 



(4.39) 



The solutions (p*+i/2i v i+i/2) f° r case ^ are summarized in the following table: 



H1-R2 


„ _ PmV m -plVl 
Pr-pl 


^i+l/2 


V i+l/2 


s > 


P« 


fj 


s < 


Pm 


Vm 


s = 


Pl+Pm 
2 


Vl+Vm 
2 



4.2.2 Computing the boundary averages from the Cauchy problem 

Liu (1979) studied the Cauchy problem for a hyperbolic system of conservation laws 
with source terms. In this subsection we apply his theory to find the boundary 
averages by solving the Cauchy problem for the PW model. 
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The homogeneous version of the PW model can be written as follows, 

u t + f(u) x = 0, (4.40) 

where u = (p,v). The Riemann problem for this system at x = has the following 
initial conditions: 

[ U h if x<0 

u(x, tj ) = { . (4.41) 

I U r , if x > 
The wave solutions to the Riemann problem consist of two basic waves with an 
intermediate state U\. By denoting U = XJ\ and U 2 = U r , we define (C^_i,C/j) as 
the i th {i = 1,2) propagation wave. Each basic wave ([/*_!, C/j) may be a shock or a 
rarefaction. These wave solutions, which have been discussed in previous subsection, 
serve as the basis for solutions to a related Cauchy problem. 

The Cauchy problem for the PW model has the following initial conditions: 

{ui(x), if x < 
(4.42) 
u r (x), if x > 

in which, 

f(ui(x)) x = s( Ul (x)), u l {Q) = U l (4.43) 

f(u r (x)) x = s(u r (x)), u r (0) = U r . (AAA) 

Here the intermediate state ui(x) is determined from 

f( Ul (x)) x = s(«i(x)), u 1 (ti) = U 1 , (4.45) 

in which U\ is the intermediate state solution to the corresponding Riemann problem. 
We denote u (x) = ui(x) and u 2 (x) = w r (x), and define (tXj_i(x), Wj(x)) as the i 4/l (i = 
1, 2) propagation wave of the Cauchy problem. 

For computational purposes, we still want to compute the boundary average of 
u(x,t) at x = 0; i.e., [/* = J u(0,t)dt. The boundary average varies when x = is 
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covered by different states or waves. Since solutions to the Cauchy problem for the 
PW model consist of three states and four types of waves, the boundary x = may 
be covered by the left, right or intermediate state, and may be crossed by 1-shock, 
2-shock, 1-rarefaction or 2-rarefaction wave. 

If the boundary x = is covered by the intermediate state, according to the 
definition of the intermediate state, the solutions of the boundary average U* is equal 
to U±. Similarly, U* = U\ or U r when the boundary is covered by the left or right 
state. 

Of the 4 types of basic waves, the 2-H and 2-R waves never cross x = since 
A 2 (£/) > 0, and 1-shock still propagates along a line described by | = cr(C/j_i, Ui). 
Therefore, these three waves don't affect the boundary averages. In the remaining 
part, we consider the adjustment to the 1-rarefaction wave. 

The 1-rarefaction wave was shown by Liu to be a perturbation on the solutions 
to the corresponding Riemann problem. Assume that (v,i-i(x) , Ui(x)) is a 1-R wave, 
Uj_i(x) and Ui(x) are separated in a region Xj_i(t) < x < Xi(t). The wave solution of 
the Cauchy problem (4.42) approaches the 1-rarefaction wave (C^_i, Ui) as t — »■ 0, 

lim^o ^r^ 1 = \i{Ui-i{x)) =^o 

x (t) ^ > 

lim^o ^r 1 = Ai(«i(s)) = 6 
lim u(x, t) — v(rj), Xi-i(t) < x < Xi(t), (4.47) 

t-»0;f=Tj 

where v(rj) G R\(Ui-i) with \i(v(r))) = rj. 

Using the coordinate of time t, we define the initial ^-characterized speed ^ as 
follows: 

^%^ = Mw&t)) (4.48) 

w(£,t) = u(x(Z,t),t) (4.49) 
and 

w(£,0) = v(0 = m (4.50) 
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Xi_i = s(£o,*) (4-51) 

u(x(£o,t),t) = Ui_i(x(^ ,t),t)=«<-i(^o,t)) ( 4 - 52 ) 

u{x{^,t),t) = u i (x(Z 1 ,t),t)=u i (x(Z 1 ,t)) (4.53) 

For a homogeneous system £ will be a constant slope for characteristics. However, £ 
is not constant for the PW model with a source term. 

With the transformation (4.48,4.49), the original system 

du(x(£,t),t) , A ^du = g ^ in w hich A(m) = V/(m) (4.54) 



dt dx 



becomes 



Since 



we obtain 



^M = Xl (w(U)) (4.55) 

dx dw , , , ,dw dx ., _„. 

aea- + (A - A % = r- (456) 



x(£,0) = 0, (4.57) 



*&o = ^=o* + o^u/ + o(* 3 ) 



<9x 1 <9 2 a; 

W t= ° t+ 2~di 

St + l^\ t = t 2 + O(?) (4.5? 



.2 , 



We calculate |rrf as follows 



at 2 

<9 2 a; 9Ai(w) du> 

St 2 " = dt 1 ~dt 

= VA lS + 0(t), (4.59) 

where we assume the difference between £|^ and A|| is small since — £wg + f(u)(- = 
for a rarefaction wave. Then, the rarefaction wave path is 

x(t) = A!t + -VAi(M)s(M)t 2 + 0(t 3 ) 
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= (--c )i+A^ 2 + O(t 3 ). (4-60) 

P 2t P 

which is a parabola. The linear term of the rarefaction wave is the 1-R wave for the 

Riemann problem. The second-order term is determined by the source term. For small 

time scales; i.e., t small, this rarefaction wave is a perturbation of the rarefaction for 

the corresponding Riemann problem. The 2-rarefaction wave functions can be derived 

similarly, although they are not necessary for computing the boundary averages. 

Recall that, by assumption, the characteristic curve of 1-rarefaction wave is 

vi = v*(p) - v*(pi). (4.61) 

Let x(t) = 0. From (4.60) we get 

v = c Q - J — 1 (4.62) 

Zrp 

I f qnr\ 

P = v- l (v-v l + v if (pi))=v~ 1 (po)-^——^- 1, (4.63) 

Wo) 2rp 

from which we can calculate U* = j Q u(0,t)dt. 

The above analysis shows that the 1-rarefaction wave fans of the Cauchy problem 
consists of parabolic curves instead of lines. When we use these parabolic rarefaction 
waves to compute U*, the numerical solution of the PW model improves. However, 
the adjustment for a 1-rarefaction wave is needed only when it crosses the boundary 
x — 0, and 1-rarefaction waves are adjusted with a lower order perturbation in a 
short time step At. Thus this improvement doesn't appear to be significant. Besides, 
we have Xi{U) = for some state U when 1-rarefaction wave cross the boundary. 
Therefore, the state U is in the unstable region for the PW model. This may be 
another reason that this adjustment doesn't yield significantly better solutions. 
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4.3 Godunov methods 

In this section, we study Godunov methods for solving the PW model (4.3). For a 
general system (4.4), the finite difference equations are 

U? +1 = U> - £(W£$) - /(Ci 1 /?)) + AW(tf), (4-64) 

in which U- U- are both averages of u(x, t) over i th cell at time (j + 1) At and j At, 
U 3 i±l L are the boundary averages calculated as shown in the preceding section, and 
s(U) is the source average over ((i — l/2)Ax, (i + 1/2) Ax) X (jAt, (j + l)At). 
When we treat the source term implicitly, the system is discretized as 

pi+ l _ p i m* 1 , /2 - m* 1 , /2 
k h 

m i+l/2 I „2*j m i-l/2 2*j 

m^ 1 — m J o* j c oPi+i/2 ~p c o/Vl/2 

ln i lri i f i+l/2 ' Pi-1/2 ' 



k h 



(4.66) 



T 

in which, p\ and ml are the cell average of p and m respectively over the ith cell, and 
p1+i/2 an< ^ ?Ti i = +i/2 are ^ ne avera g es of p and m respectively on the cell boundary x i+ i/ 2 
in the time interval (tj,tj+i). In (4.66), the source term s(U) is treated implicitly. 
From (4.65,4.66), we can write the evolution equations for the PW model as 
k 

pl +1 =pI~ r( m *ii/ 2 - m *-i/ 2 ) ( 4 - 67 ) 



v\""i+l/2 "*i-l/2> 

n i+l/2 

Pi+1/2 Pi+1/2 



*7 / *1 

1 , „■ k™' /2 „ . (m. J 



j+l l r j K,"*i+l/2 2 *i \""i-l/2J 2 *j I /,, «o\ 

m * = ?TTI^ m * ~ /^~^ + c o^ii/2-— 7T c oP/i/ 2 ] ( 4 - 68 ) 



+-/;(p! +1 )} 

r 
In a first-order Godunov method we use the cell averages p\ — pi, mi = ml and 
= pl +1 ,m r = m 3 i+1 as left/right states as the initial condition for the Riemann 
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problem 

Ui, if x — Xj+i/2 < 
u i+1 /2(x,tj) = { ' . (4.69) 

U r , if x - Xj+i/2 > 
4.3.1 The Second-order Godunov Method 

In this subsection we introduce a second-order Godunov method for the PW model 
in the process similar to that in section 3.3. 

We begin by writing the PW model in the linearized form 

u t + A(u)u x = s(u), (4.70) 

where 

u = |' PiX ' t] ] (4.71) 

m(x, t) 



and 



AW) = 1 ° * I. (1.72) 



m 2 I 2 2m 
2 + c 



The eigenvalues and eigenvectors of A(u) are 

A!(w) = * - Co, rx(«) = [1, AJ* = [1, f - c ]' 
A 2 (w) = s + Co , ri ( u ) = [i, A 2 ]* = [1, s + Co ]* 

We diagonalize A(u) by 



where the transformation matrix T(u) is 



(4.73) 



r^uji^^u) = | | =A(u), (4.74) 

A 2 («) 



T{u) = \ \. (4.75) 



in in I 

7 - c 7 + co 
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Letting W = T^ 1 (u)u, the PW model under the transformation becomes 

W t + K{u)W x = T-\u)s(u). (4.76) 

Therefore, the PW model is transformed into two separated scalar equations in W = 
(wi,W2). For the scalar equation in Wi (i = 1,2), we introduce an interpolation for 
Wi(x,t) which yields a second-order Godunov method for solving this equation. In 
the remaining part of this subsection, we first introduce an interpolation for Wi(x,t) 
and then apply inverse transformation on them in order to develop a second-order 
method for the whole system. 

For a scalar equation wt + X(w)w x = 0, in a first-order Godunov method we use 
a step function ui(x,tj) to interpolate the data at time tj, 

w I (x,t j ) = wf, if Xj_i/ 2 < x < Xj+i/2, (4.77) 

and solve the Riemann problem with the following initial conditions: 



w i+i/2 = < 

i,R i 

W i-l/2 = W i 



(4.78) 



In a second-order Godunov method, we interpolate the data at time tj with a 
piecewise linear function, 

w T (x, tj) = wf + ^ X ~ l > A VL wl if x 4 _ 1/2 < x < x i+1/2 . (4.79) 

With a half-step prediction , the Riemann problem has the initial conditions of 

<ll} L = ^ + |(i-ak)!)a^X (480) 

!££$•* = ^-i(i + AK)f)A^v' 

where A VL w J i is the van Leer slope (all superscripts j are suppressed) 

S i -mm(2\w i+1 -Wi\,2\w i -w i . 1 \,\\w i+1 -w i _ 1 \), (p>0 
A vn Wi = \ 

0, otherwise 
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Si = sign(w i+1 - wt-x) (4.81) 

(p = (w i+ i - Wi) ■ (wi - Wi-i) (4.82) 

We apply the above procedure twice to compute the half-step values W? +1 L ' 
and W-_ 1/2 ' . Then we apply the inverse transformation on these half-step values to 
obtain U 3 i+l j 2 ' and U^_ 1 L ' : 



rrj + l^L _ ( j. j + l/2,L 

U i+l/2 ~ A \ U i ) VV i+l/2 

T j+1/2,R _ , js j+l/2,R 



(4.83) 



With this new interpolation of p and m, we can solve the Riemann problem or 
the Cauchy problem on a cell boundary. For the inhomogeneous system, this second- 
order method has a convergence rate of 2. However, it is different for the PW model 
with a relaxation term. 

4.3.2 Some other Godunov-type variant methods 

In this subsection, we review other variants of Godunov method for (4.3) with a 
source term. 

The first variant was suggested by Pember (1993a, 1993b). He treated the source 

term as s(U) = U^i-i/^ + s ( U !+i/2^ where both ^-V? and U i+W are the 
boundary averages solved in the Riemann problem. 

The second variant is the fractional step splitting method, in which each time 

step At is split in-to three steps. In the first and third fractional steps, a first-order 

implicit method is used to integrate 





(4.84) 
for time steps of At/2. In the second step, we solve the corresponding homogeneous 
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system of (4.3), i.e., 

P ] + [ 2 m 1=0 (4.85) 

for a time step of At. 

A third variant is the quasi-steady wave-propagation algorithm suggested by 
LeVeque (1998a, 1998b). This method introduces a new discontinuity in the center of 
each grid, i.e., 

l(U~ + U?) = U u (4.86) 

and 

f(Ut)-f(UT) = 'WAX, (4.87) 

i Pt \ 
where U i — \ . Then we get 

mf j 

mf - m~ \ ( 

= Ax (U 

so mf = m~ = rrii. Setting 

pf = p t + 5, p~ = Pi -5 (4.89) 

we find 

2c 2 5 3 - Kb 2 - (2m 2 + 2c 2 p 2 )5 + Kp 2 = 0, (4.90) 

in which K = ■'*^ pl '~ mi Ax. 

T 

Given the solution to (4.90), LeVeque solved the Riemann problem of the homo- 
geneous system with the initial conditions 

, Uf if x — Xi+i/o < 
ik+iftfatj) = { l +/ (4.91) 

U i if x - x i+1/2 > 
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Then in the evolution equations, the source term is not considered. 

In LeVeque's method, the data at time tj are interpolated with a solution to the 
standing wave for the PW model in each cell. This method may be used together 
with the Cauchy problem we discussed before. 

4.4 Numerical Solutions to the PW Model 

In this section, we use the model parameters given in (Kerner, 1994), i.e., cq = 
2.48445//r, u*(p) = V(p) = 5.0461[(1 + exp{[p-0.25]/0.06})- 1 -3.72 x Kr 6 ]Z/r and 
L = 800/. Here / is the unit of length, r is the relaxation coefficient, and the section 
of the roadway is from 0/ to 800/. We set / = 10 (m) r = 10 (sec) and ph = 0.172. 
The equilibrium functions i>*(p) and /*(p) are given in Figure 4.1. In the figure, 
Pel, Pci are two critical densities and the region p c i < p < p C 2 is the unstable region. 
According to (4.9), p ci and p c2 satisfy the equation: 

pv'M+co = 0. (4.92) 

From this equation we get p c \ = 0.173,p C 2 = 0.396. 

For numerical computation purpose, we use two initial conditions: 

p(x,0) = p fe + 0.02sin(27Tx/800//) (4.93) 

v(x,0) = v,{p h )-0.02cos{27ix/800/l), (4.94) 

which is a global perturbation, and 

p h + 5 when x G [37.5/, 48.4/] 

p(x, 0) = I p h - 5/3 when x G [50.0/, 82.8/] (4.95) 

Ph otherwise 

v 

v(x,0) = v*(p(x,0)), (4.96) 

which is a local perturbation. Setting S = 0.02 and ph = 0.16, the two initial condi- 
tions are given in Figure 4.2. The global perturbation is not in equilibrium, however 
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the local perturbation is in equilibrium. 

4.4.1 Stability test 

The PW model is unstable in some regions, which is different from Zhang's model. In 
this subsection we test the stability property of the PW model at time T = 200. Or. 
We use a first-order Godunov's method to compute the relative differences, defined 
in (3.73), between solutions with N = 512 grids and those with 1024 grids, and 
the difference between solutions with 1024 grids and those with 2048 grids. These 
differences are drawn as curves labeled as 1024 — 512 and 2048 — 1024 respectively in 
the following figures. Since the solutions (p, v) are close to the equilibrium state, only 
the differences of p are given. Given the convergent Godunov's method, the difference 
caused by different grid numbers decreases if the PW model is stable. We test the 
PW model with different initial conditions and different boundary conditions, which 
are in the stable or unstable region for the PW model, and the results are the same 
as predicted. 

Setting ph = 0.16 or p^ = 0.17 for the initial conditions (4.93,4.94), we solve the 
PW model with 512, 1024 and 2048 grids with periodic boundary conditions. The 
differences are shown in Figure 4.3. In the figure, the difference decreases when 
we increase the number of grids when ph = 0.16; but the difference increases when 
Ph = 0.17. The figure proves that the PW model is stable for ph = 0.16, and unstable 
for p h = 0.17. 

Using the same initial conditions as in last example, we solve the PW model with 
Neumann boundary conditions. The relative differences are given in Figure 4.4. In 
this case we get the same results as in last example. 

Next we use the initial conditions (4.95,4.96) with 8 = 0.02. We solve the PW 
model with periodic boundary conditions. According to the differences shown in 
Figure 4.5, we conclude that the PW model also stable when ph = 0.16 and unsta- 
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Figure 4.1: One selection of the equilibrium velocity and flow rate 
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Figure 4.2: Two initial conditions: global non-equilibrium perturbation v.s. local equilibrium per- 
turbation 
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Figure 4.3: Stability and instability for (4.93,4.94) with Periodic boundary condition 
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Figure 4.4: Stability and instability for (4.93,4.94) with Neumann boundary condition 
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Figure 4.5: Stability and instability for (4.95,4.96) with Periodic boundary condition 
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R1-R2: evolution of p (density) 
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R1-R2: evolution of v(I/t) 
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Figure 4.6: The R1-R2 wave solutions to the Ricmann problem 

ble when ph = 0.17. We don't test the PW model for the local perturbation with 
Neumann boundary conditions, since we expect the same stability property. 

According to the numerical tests, the system is stable when maxp(x, 0) < 0.18 
and unstable when maxp(x, 0) > 0.19, which is consistent with the theoretical pre- 
diction. 



4.4.2 The Riemann problem and steady-state solutions 

In this subsections, we solve the PW model with four well-selected initial conditions 
so that we observe second-order waves. The Neumann boundary conditions are used 
and the number of grids is set as N = 1024. We also change the relaxation term 
to , " m . The relaxation time lOOOr is long enough for us to watch the second-order 

IOOOt ° ° 

waves at T n = lOOr, before all 2-waves relax to 1-waves. 
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R1-H2: evolution of p (density) 
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Figure 4.7: The R1-H2 wave solutions to the Riemann problem 
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Figure 4.8: The H1-H2 wave solutions to the Riemann problem 
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H1-R2: evolution of p (density) 
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Figure 4.9: The H1-R2 wave solutions to the Ricmann problem 



1. We use the following jump initial conditions 



p(x,0) 

v(x, 0) 



0.16 

vi = t>*(0.16) when x <= 400/ 

v r = ^(0.16) + 0.21/t when x > 400/ 

The solutions to the Riemann problem are a R1-R2 wave, given in Figure 4.6. 



(4.97) 
(4.98) 



2. We use the following jump initial conditions 



p(x,0) 
v(x, 0) 



pi = 0.16 when x <= 400/ 

p r = 0.16 -0.02 when x > 400/ 
^(0.16). 



(4.99) 
(4.100) 



The wave solutions to the Riemann problem of the PW model are a R1-H2 wave, 
given in Figure 4.7. 
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3. We use the following jump initial conditions 

p(x,0) = 0.16 (4.101) 

{ vi = u»(0.16) when x <= 400/ 

v(x,0) = { . (4.102) 

[v r = u,(0.16) - 0.2//r when x > 400/ 
The wave solutions are a H1-H2 wave, given in Figure 4.8. 

4. We use the following jump initial conditions 

| pi = 0.16 when x <= 400/ 

p(x,0) = { (4.103) 

{ p r = 0.16 + 0.02 when x > 400/ 
v(x,0) = u,(0.16). (4.104) 

The wave solutions are a H1-R2 wave, given in Figure 4.9. 

The solutions above show four different type of second-order waves, which consist 
of two basic waves. However, 2- waves relax to 1-waves if the relaxation time is shorter 
or the observing time is longer. A 2-shock wave relaxes to a 1-rarefaction wave; and 
a 2-rarefaction wave relaxes to a 1-shock wave, due to the effect of the source term. 
In next part, we show how the 2-waves relax to 1-waves and how the free regions and 
cluster regions form. We set relaxation time as lOr, and observe the solutions until 
T = lOOr. 

1. With the initial conditions (4.97, 4.98), we get the solutions shown in Figure 

4.10. At around 5r a downstream 1-shock forms when the traffic conditions 
relax to the equilibrium state, i.e., v = i>*(0.16). After that traffic flow forms a 
free region with lower density and higher travel speed. The free region travels 
in the speed of A*. However the free region will disappear as the Rl-Hl wave 
propagates, and finally the traffic flow will become uniform. 

2. With the initial conditions (4.99, 4.100), we get the solutions shown in Figure 

4.11. A new 1-rarefaction wave forms when the H2-wave disappears. As these 
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two rarefaction waves propagate, the traffic conditions will become uniform. 

3. With the initial conditions (4.101, 4.102), we get the solutions shown in Figure 

4.12. At around 5r, the 2-shock wave disappears and a 1-rarefaction wave forms. 
After that a cluster, with higher density and lower travel speed, travels with at 
the speed of A*. However, the traffic conditions will become uniform as the 
Hl-Rl wave propagates. 

4. With the initial conditions (4.103, 4.104), we get the solutions shown in Figure 

4.13. As long as the 2-rarefaction disappears, a new 1-shock is formed. After 
that, both shock waves travel in the speed of A*(0.16) = 2.121 /t. 

In the remaining part of this subsection we consider the steady-state solutions of 
the system after a long time T = 1600r with the relaxation time r . 

1. Using the initial conditions (4.93,4.94) with ph = 0.16 and periodic boundary 
conditions, we get the solutions shown in Figure 4.14. The figure shows that 
an upstream shock wave and a downstream rarefaction wave form. After that, 
the traffic conditions become more and more uniform. 

2. Using the initial conditions (4.95,4.96) with p^ = 0.16 and periodic boundary 
conditions, we have the solutions shown in Figure 4.15. The solutions show that 
the traffic conditions get more and more uniform, quicker than the case shown 
in Figure 4.14. 

4.4.3 General solutions and convergence rates 

In this subsection, different numerical methods for the PW model are discussed with 
the initial conditions (4.93,4.94) (p^ = 0.16). Using Neumann boundary conditions, 
we solve the PW model until T = 400r. With the number of grids as 64, 128, 256, 512 
or 1024, we carry out 5 separate computations. For 1024 grids, we plot the contour 
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Evolution of p (density) when relaxation time = 5 x 
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Figure 4.10: Formation of a free region 
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Evolution of v (I/ t) when relaxation time = 5 z 
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Figure 4.11: Formation of two 1-rarefaction waves 
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Evolution of p (density) when relaxation time = 5 x 
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Evolution of v (l/x) when relaxation time = 5 x 
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Figure 4.12: Formation of a cluster 



100 

80 

. 60 

: 40 

20 




Evolution of p (density) when relaxation time = 5 x 




100 200 300 400 500 600 700 

x/l 



Evolution of v (l/x) when relaxation time = 5 t 




4.1 

4.05 

4 

3.95 

3.9 

3.85 



Figure 4.13: Formation of two 1-shock waves 
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p 


128-64 


Rate 


256-128 


Rate 


512-256 


Rate 


1024-512 


L 1 


7.36e-04 


9.72e-01 


3.75e-04 


9.95e-01 


1.88e-04 


9.69e-01 


9.62e-05 


L 2 


9.38e-04 


6.65e-01 


5.92e-04 


6.81e-01 


3.69e-04 


6.74e-01 


2.31e-04 


L°° 


3.07e-03 


4.02e-02 


2.99e-03 


1.79e-01 


2.64e-03 


3.27e-01 


2.10e-03 


V 


128-64 


Rate 


256-128 


Rate 


512-256 


Rate 


1024-512 


L 1 


9.68e-03 


9.67e-01 


4.95e-03 


9.99e-01 


2.48e-03 


9.59e-01 


1.27e-03 


L 2 


1.27e-02 


6.43e-01 


8.14e-03 


6.63e-01 


5.14e-03 


6.56e-01 


3.27e-03 


L°° 


4.41e-02 


4.81e-02 


4.27e-02 


1.76e-01 


3.78e-02 


3.15e-01 


3.04e-02 



Table 4.1: Convergence rates for the first-order Godunov's method 

of p and v on the x — t phase plane as well as their solutions at selected times. We 
also compute the convergence rates based on solutions with different number of grids. 
The convergence rate is defined in (3.73 - 3.74). 

The first-order Godunov's method gives the solutions shown in Figure 4.16 and 
Figure 4.17. The relative errors and convergence rates are given by Table 4.1. From 
the table, we can see that for L 1 norm the method is almost of first order, but for L 2 
or L°° norms, the rate of convergence is lower. 

The quasi-steady wave-propagation algorithm by LeVeque (1998a, 1998b) gives 
the solutions shown in Figure 4.18 and Figure 4.19 for 1024 grids. The relative 
errors and convergence rates are given in Table 4.2. This scheme de-estimate the 
effects of the source term, since the convergence rates of p and v are totally different. 

The fractional step splitting method gives solutions of the PW model shown in 
Figure 4.20 and Figure 4.21 for 1024 grids. The relative errors and convergence 
rates are given in Table 4.3. We find that the fractional step splitting method is not 
so stable as the first-order method, since convergence rates have big oscillations. 
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Evolution of p (density) 
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Figure 4.14: Solution for (4.93,4.94) till 1600t 
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Figure 4.15: Solution for (4.95,4.96) till 1000t 
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Evolution of p (density) with Neumann BC 
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Figure 4.16: Solutions by a first-order Godunov's method for 1024 grids 
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Figure 4.17: Solutions from Figure 4.16 at selected times 
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Evolution of p (density) with the quasi-steady wave-propogation algorithm 
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Figure 4.18: Solutions by LeVeque's method for 1024 grids 



solutions of p at different times 




Figure 4.19: Solutions from Figure 4.18 at selected times 



CHAPTER 4. 



The PW Model and Its Numerical Solutions 



73 



p 


256-128 


Rate 


512-256 


Rate 


1024-512 


L 1 


2.09e-02 


5.10e-01 


1.47e-02 


1.16e+00 


6.54e-03 


L 2 


2.28e-02 


4.88e-01 


1.63e-02 


1.14e+00 


7.40e-03 


L°° 


3.49e-02 


5.81e-01 


2.33e-02 


8.86e-01 


1.26e-02 


V 


256-128 


Rate 


512-256 


Rate 


1024-512 


L 1 


3.42e-02 


-1.61e+00 


1.05e-01 


6.18e-02 


1.00e-01 


L 2 


3.91e-02 


-1.58e+00 


1.17e-01 


5.78e-02 


1.12e-01 


L°° 


6.96e-02 


-1.30e+00 


1.71e-01 


-1.51e-01 


1.90e-01 



Table 4.2: Convergence rates for the quasi-steady wave-propagation algorithm 



Evolution of p (density) with fractional step splitting method 
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Figure 4.20: Solutions by fractional splitting method for 1024 grids 
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solutions of p at different times 




Figure 4.21: Solutions from Figure 4.20 at different times 



p 


128-64 


Rate 


256-128 


Rate 


512-256 


Rate 


1024-512 


L 1 


8.27e-04 


9.49e-01 


4.28e-04 


-5.45e-01 


6.25e-04 


1.38e+00 


2.41e-04 


L 2 


1.00e-03 


7.35e-01 


6.01e-04 


-2.72e-01 


7.26e-04 


5.81e-01 


4.85e-04 


L°° 


2.83e-03 


1.95e-01 


2.47e-03 


-4.07e-01 


3.28e-03 


-4.48e-01 


4.47e-03 


V 


128-64 


Rate 


256-128 


Rate 


512-256 


Rate 


1024-512 


L 1 


1.08e-02 


9.46e-01 


5.58e-03 


-5.01e-01 


7.90e-03 


1.35e+00 


3.11e-03 


L 2 


1.33e-02 


7.14e-01 


8.13e-03 


-1.92e-01 


9.28e-03 


4.76e-01 


6.67e-03 


L°° 


3.86e-02 


1.19e-01 


3.55e-02 


-2.77e-01 


4.30e-02 


-5.32e-01 


6.22e-02 



Table 4.3: Convergence rates for the fractional step splitting method 
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Evolution of p (density) with another treatment of source term 
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Figure 4.22: Solutions by Pember's method for 1024 grids 



solutions of p at different times 




Figure 4.23: Solutions from Figure 4.22 at different times 
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p 


128-64 


Rate 


256-128 


Rate 


512-256 


Rate 


1024-512 


L 1 


8.31e-04 


7.69e-01 


4.88e-04 


9.73e-01 


2.49e-04 


1.02e+00 


1.22e-04 


L 2 


1.07e-03 


4.67e-01 


7.76e-04 


6.39e-01 


4.98e-04 


6.82e-01 


3.10e-04 


L°° 


2.94e-03 


-4.29e-02 


3.03e-03 


9.29e-02 


2.84e-03 


2.85e-01 


2.33e-03 


V 


128-64 


Rate 


256-128 


Rate 


512-256 


Rate 


1024-512 


L 1 


1.10e-02 


7.60e-01 


6.50e-03 


9.73e-01 


3.31e-03 


1.02e+00 


1.64e-03 


L 2 


1.46e-02 


4.52e-01 


1.06e-02 


6.33e-01 


6.87e-03 


6.74e-01 


4.30e-03 


L°° 


4.09e-02 


-4.25e-02 


4.21e-02 


9.82e-02 


3.94e-02 


2.78e-01 


3.24e-02 



Table 4.4: Convergence rates for Pember's method 

Pember's method (1993a, 1993b) give the solutions shown in Figure 4.22 and 
Figure 4.23 for 1024 grids. The relative errors and convergence rates are given in 
Table 4.4. 

In the remaining part of this subsection we consider the second-order Godunov's 
method. We use ph = 0.15 for initial conditions (4.93,4.94), since the second-order 
method is not stable for ph = 0.16. For 1024 grids, we get the solutions shown in 
Figure 4.24 and Figure 4.25. The relative errors and convergence rates are given in 
Table 4.5. 

Comparing the convergence rates with those for first-order Godunov's method, 
we see no significant improvement. This is different from the case for Zhang's model. 
This is a special property of the PW model. 

4.4.4 Unstable solutions of the PW model 

In this subsection we check the unstable solutions of the PW model. We use the 
first-order Godunov's method with periodic boundary conditions, and we observe the 
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Evolution of p (density) with second order method 






0.16b 




0.16 
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Figure 4.24: Solutions by the second-order Godunov's method for 1024 grids 



solutions of p at different times 




solutions of v corresponding to p 
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Figure 4.25: Solutions from Figure 4.24 at selected times 
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p 


128-64 


Rate 


256-128 


Rate 


512-256 


Rate 


1024-512 


L 1 


5.34e-04 


1.29e+00 


2.18e-04 


1.19e+00 


9.55e-05 


1.03e+00 


4.69e-05 


L 2 


5.34e-04 


1.29e+00 


2.18e-04 


1.19e+00 


9.55e-05 


1.03e+00 


4.69e-05 


L°° 


5.34e-04 


1.29e+00 


2.18e-04 


1.17e+00 


9.65e-05 


1.01e+00 


4.78e-05 


V 


128-64 


Rate 


256-128 


Rate 


512-256 


Rate 


1024-512 


L 1 


6.02e-03 


1.30e+00 


2.45e-03 


1.19c+00 


1.07e-03 


1.03e+00 


5.26e-04 


L 2 


6.02e-03 


1.30e+00 


2.45e-03 


1.19e+00 


1.07e-03 


1.03e+00 


5.26e-04 


L°° 


6.02e-03 


1.30e+00 


2.45e-03 


1.18e+00 


1.08e-03 


1.01e+00 


5.36e-04 



Table 4.5: Convergence rate for second-order method 

solutions at To = 200r for different number of grids. 

We use the initial conditions (4.93,4.94) with ph = 0.175, which is in the unstable 
region of the PW model. Solutions of the PW model are listed as the following for 
different number of grids: 

1. For 512 grids, solutions are given in Figure 4.26 and Figure 4.27. 

2. For 1024 grids, solutions are given in Figure 4.28 and Figure 4.29. 

3. For 2048 grids, solutions are given in Figure 4.30 and Figure 4.31. 

The figures show that the number and position of spikes are different for differ- 
ent number of grids. This difference is caused by different approximations used by 
Godunov's method for different number of grids, since the PW model is unstable in 
the region where the initial conditions are. 
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Unstable evolution of p (density) when N=512 




100 200 300 400 500 600 700 



Unstable evolution of v (I/ -r) when N=51 2 




100 200 300 400 500 600 700 



Figure 4.26: Solutions for 512 grids with initial conditions (4.93,4.94) 



solutions of p at different times 




Figure 4.27: Solutions from Figure 4.26 at selected times 
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Unstable evolution of p (density) when N=1024 




4 
3.5 

3 

2.5 

2 

1.5 

1 



100 200 300 400 500 600 700 



Figure 4.28: Solutions for 1024 grids with initial conditions (4.93,4.94) 



solutions of p at different times 




Figure 4.29: Solutions from Figure 4.28 at selected times 
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Unstable evolution of p (density) when N=2048 
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Figure 4.30: Solutions for 2048 grids with initial conditions (4.93,4.94) 



solutions of p at different times 










solutions of v corresponding to p 
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Figure 4.31: Solutions from Figure 4.30 at selected times 



Chapter 5 

The Inhomogeneous LWR Model 
and Its Numerical Solutions 



5.1 Introduction 

The LWR model (Lighthill and Whitham, 1955;Richards, 1956) was introduced based 
on the conservation of traffic flow, and is written as: 

£> + i' = »■ '"J 

The LWR model assumes that traffic flow is in equilibrium, or equivalently that 
travel speed v (unit: mph) is defined as a function of traffic density p (unit: vpm) at 
any location x: 

v = v m (x,p). (5.2) 

The function of flow rate / = pv*(x,p) (unit: vph) is called a fundamental diagram. 
For typical equilibrium traffic flows, travel speed is a decreasing function; i.e., v* p < 0, 
and traffic flow rate is a concave function; i.e., f pp < 0. 

The homogeneous LWR model is for modeling a homogeneous road, where travel 
speed t>* is uniform with respect to location x. The homogeneous LWR model can 

82 
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thus be written as 

Pt + f(p)x = 0. (5.3) 

The homogeneous LWR model is a scalar conservation law, and there have been 
many methods to compute its entropy solutions (Lebacque, 1995). It's well-known 
that the entropy solutions exist and are unique under the so-called "Lax entropy 
condition". For computation of these entropy solutions, the Godunov method is 
most often used. 

For a road with inhomogeneity, such as variable number of lanes, curvature and 
slope, we can formulate the inhomogeneous LWR model as 

Pt + f(a,p) x = 0, (5.4) 

where a = a(x) is a time-invariant variable. The inhomogeneity factor a(x) gives a 
profile of a piece of roadway, for example a(x) can be the number of lanes at location 
x. 

The inhomogeneous LWR model has been studied by Daganzo (1994) and Lebacque 
(1995) and Daganzo (1994). Both of these authors suggested solutions to the inhomo- 
geneous LWR model, and their solutions are consitent. However, these studies only 
presented empirical solutions without rigorous proof. 

The difficulty of dealing with the inhomogeneous LWR model is due to the ex- 
tra variable a(x). Here, by introducing a(x) as an additional conservation law, we 
consider the inhomogeneous LWR model as a resonant nonlinear system, which has 
been discussed in (Isaacson & Temple, 1992; Lin et al., 1995). In this chapter, we 
follow the procedures provided in those researches to study the inhomogeneous LWR 
model. 

This chapter is organized as follows. In section 2, we formulate the inhomogeneous 
traffic model as a resonant nonlinear system, and its properties are discussed. In 
section 3, we solve Riemann problem for the inhomogeneous LWR model. In section 
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4 we present the numerical methods for this model. We conclude the discussions of 
this chapter in section 5. 



5.2 The Properties of the inhomogeneous LWR model 

To apply the results for hyperbolic systems of conservation laws, we express the 
inhomogeneity factor a(x) as an additional conservation law, i.e., a t = 0. Hence, we 
can write the inhomogeneous LWR model as 



U t + F(U) x = 0, 



(5.5) 



where U = (a,p),F(U) = (0, f(a, p)),x E R, t > 0. In this chapter we consider one 
type of inhomogeneity - variable number of lanes. The fundamental diagram is thus 
written as f(a,p) = pf*(-), given that all the lanes are of the same condition. 
The inhomogeneous LWR model (5.5) can be linearized as 



U t + DF(U)U X = 0, 
where the differential DF(U) of the flux vector F(U) is 



(5.6) 



DF 



The two eigenvalues for DF are 



-Pl v '(£. 



W(S) «.(£) + MS) 



(5.7) 



A = \ l=v ^) + H v '^) 



P ,,Ps 

v * 

a a' 



(5. 



The corresponding right eigenvectors are 



Ro 



«.(£) + MS) 



EYv'(P) 



Ri 





1 



and the left eigenvector of df /dp as L = 1. The system (5.5) is a non-strictly 
hyperbolic system, since Ai may be equal to Aq = 0. 



CHAPTER 5. The Inhomogeneous LWR Model and Its Numerical Solutions 



85 



We consider a traffic state [/* = (a*,p*) as critical if 



(5.9) 



i.e., at critical states, the two wave speeds are the same and system (5.5) is singular. 
For a critical traffic state [/* we also have 

d 



dp 



Xi(K 



fop < o, 



and 



Ijw 



:^) 2 v'jP) 



I £7, 



-v*(-)\ Uf >0. 



(5.10) 



(5.H) 



a a a a 

A consequence of properties (5.10) and (5.11) is that the linearized system (5.6) 
at (7* has the following normal form 



<5a 


+ 





o~a 


. ^. 


t 


1 


_^_ 



0. 



(5.12) 



The system (5.12) has the solution 5p(x,t) = 5a'(x)t + c, and the solution goes to 
infinity as t goes to infinity. Therefore (5.12) is a linear resonant system, and the 
original inhomogeneous LWR model (5.5) is a nonlinear resonant system. 

For (5.5), the smooth curve T in [/-space formed by all critical states [/* are 
named a transition curve. Therefore T is defined as 



r 



{I7|A 1 (17)=0}. 



Since Ai(J7) = v*(f ) + f<(f ), we obtain 

T = < (a, p)\— = a, where a uniquely solves v*(a) + av'^a) = > ; (5.13) 

i.e., the transition curve for (5.5) is a straight line passing through the origin in 
(7-space. In (5.13), a is unique since f(a,p) is concave in p. 

The entropy solutions to a nonlinear resonant system are different from those to 
a strict hyperbolic system of conservation laws. Isaacson & Temple (1992) proved 
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0-wave 



1-wave: a=a* 



500 



1000 



1500r 



Figure 5.1: Integral curves 



that solutions to the Riemann problem for system (5.5) exist and are unique with 
the conditions (5.9)-(5.11). Lin et al. (1995) also presented solutions to a scalar 
nonlinear resonant system, which is similar to our system (5.5) except that / is 
convex. In the next section we apply those results to solve the Riemann problem for 
the inhomogeneous LWR model. 
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5.3 Solutions to the Riemann problem 

In this section we study the wave solutions to the Riemann problem for (5.5) with 
the following jump initial conditions 

[ U L if x < 
U(x,t = 0) = { , (5.14) 

[ U R iix>0 

where the initial values of Ul, Ur are constant. For computational purpose, we are 
interested in the average flux at the boundary x = over a time interval At, which 
is denoted by f^; i.e., we want to find 

1 /-A* 



fS = At/ o f(U(x = 0,t))dt. (5.15) 

The inhomogeneous LWR model (5.5) has two families of basic wave solutions 
associated to the two eigenvalues. The solutions whose wave speed is Ao are in the 0- 
family, and the waves are called 0-waves. Similarly the solutions whose wave speed is 
Ai are in the 1-family, and the waves are called 1-waves. In [/-space, the wave curves 
for (5.5) are the integral curves of the right eigenvectors Ro and R4. Hence the 0-wave 
curves are given by f(U)— const, and the 1-wave curves are given by a = a, where a is 
constant. The 0-wave is also called a standing wave since its wave speed is always 0. 
The 1-wave solutions are determined by the solutions of the scalar conservation law 
Pt + f{a,p)x = 0. A 0-wave curve, a 1-wave curve passing a critical state [/* and the 
transition curve F are shown in Figure 5.1, where a is set as the vertical axis and p 
is set as the horizontal axis. 

As shown in Figure 5.1, the 0-wave curve is convex, and the 1-wave curve is 
tangent to the 0-wave curve at the critical state [/*. The transition curve Y intersects 
the 0-wave and 1-wave curves transversely at [/*, and there is only one critical state 
on one 0-wave or 1-wave curve. For any point U, there is only one 0-wave curve and 
only one 1-wave curve passing it. In Figure 5.1, the states left to the transition 
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a 3- 



Figure 5.2: The Riemann problem for Ul left of T 

curve are undercritical since p/a < a; and the states right to the transition curve are 
overcritical since p/a > a. 

The wave solutions to the Riemann problem for (5.5) are combinations of basic 
0-waves and 1-waves. If we order the waves with respect to space x at any time t, all 
the waves must satisfy Lax's entropy condition; i.e., the waves from left (upstream) 
to right (downstream) should increase their wave speeds so that they don't cross 
each other. This condition is imposed on all hyperbolic systems of conservation laws. 
For a nonlinear resonant system (5.5), an additional condition has to be imposed; 
i.e., as long as the standing wave is not interrupted by a shock- wave, its associated 
density must vary continuously. To guarantee this, the following entropy condition is 
required: 



The standing wave can NOT cross the transition curve T. 



(5.16) 



With the two entropy conditions, the solutions to the inhomogeneous LWR model 
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Figure 5.3: The Riemann problem for Ul right of T 

exist and are unique. The wave solutions for undercritical left state Ul is shown in 
Figure 5.2, and those for overcritical left state Ul is shown in Figure 5.3. 

In the remaining part of this section, we discuss wave solutions to the Riemann 
problem for (5.5), present the formula for the boundary flux /q related to each type of 
solutions, summarize our results and compare them with those existing in literature. 



5.3.1 Solutions of the boundary fluxes 

When Ul = (cll,Pl) is undercritical; i.e., Pl/ul < a , where a is defined in (5.13), 
we denote the special critical point on standing wave passing Ul as [/*. Thus, as 
shown in Figure 5.2, the ?7-space is partitioned into three regions by DU*, OU* 
and C/*C, where DU* = {(a, p)\a = a*, p < p*}, OU* = F fl {0 < p < p*} and 
C/*C = {(a,p)|/(a,p) = J'{Ul),P > p*}- Related to different positions of the right 
state Ur in the [/-space, the Riemann problem for (5.5) with initial conditions (5.14) 
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Figure 5.4: One example for wave solutions of type 1 for (5.5) with initial conditions (5.14) 
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Figure 5.5: One example for wave solutions of type 2 for (5.5) with initial conditions (5.14) 
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Figure 5.6: One example for wave solutions of type 3 for (5.5) with initial conditions (5.14) 
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Figure 5.7: One example for wave solutions of type 4 for (5.5) with initial conditions (5.14) 



CHAPTER 5. The Inhomogeneous LWR Model and Its Numerical Solutions 92 

has the following six types of wave solutions. After discussion for each type of solutions 
we provide formula for calculating the associated boundary flux /q. 

Type 1 When U R is in region ABUiJJ^DA shown in Figure 5.2; i.e., 

f(U R ) < /(£/*) = f(U L ), p R /a R < a and a R > a*, (5.17) 

wave solutions to the Riemann problem are of type 1. These solutions consist of 
two basic waves with an intermediate state U\ = {a R , Pi\f( aR , Pl )=f(u*)=f(u L ))- Of 
these two waves, the left (Ul, Ui) is a standing wave, and the right (U±, U R ) is a 
rarefaction wave with characteristic velocity X\(a, p) > 0. 

From Figure 5.2, we can see that the Riemann problem may admit this type of 
solutions when a^ > a R or a^ < a R ; i.e., when the road merges or diverges at 
x = 0. Here we present an example of this type of solutions in Figure 5.4, where 
the roadway merges at x = 0; i.e., az, > a R - For the case when the roadway 
diverges at x — 0, we can find similar solutions. 

From Figure 5.4, we obtain the boundary flux /q = fiU^) = f(U*) for wave 
solutions of type 1. 

Type 2 When U R is in region BUlU^CB shown in Figure 5.2; i.e., 

f(U R ) > f(U.) = f(U L ), (5.18) 

wave solutions to the Riemann problem are of type 2. These solutions consist of 
two basic waves with an intermediate state U\ = (<1r, Pi\f( aR , pi )=f(u*)=f(u L ))- Of 
these two waves, the left (Ul, Ui) is a standing wave, and the right (Ui, U R ) is a 
shock wave with positive speed a = — R '~^ *> > 0. 

From Figure 5.2, we can see that the Riemann problem may admit this type of 
solutions when the downstream traffic condition U R is undercritical or overcrit- 
ical, or the roadway merges or diverges at x — 0. Here we present an example 
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of this type of solutions in Figure 5.5, where the downstream traffic condition 
overcritical and the roadway merges at x — 0. For other situations as long as 
(5.18) is satisfied, we can find similar solutions. 

From Figure 5.5, we obtain the boundary flux /q = f{Ui) = /(£/*) for wave 
solutions of type 2. Here we have the same formula as that for wave solutions of 
type 1. 

Type 3 When U R is in region OUJJO shown in Figure 5.2; i.e., 

f(U R ) < /([/*) = f(U L ), p R /a R > a, (5.19) 

wave solutions to the Riemann problem are of type 3. These solutions consist 
of two basic waves with an intermediate state U\ = (a>L, Pi\f(a L , Pl )=f(u R )) ■ Of 
these two waves, the left one (Ul, Ui) is a shock wave with negative speed o = 
•^ 1 _ ' L ' < 0, and the right one (U\, U R ) is a standing wave. 

From Figure 5.2, we can see that the Riemann problem may admit this type 
of solutions when the roadway is merges or diverges at x — 0. Here we present 
an example of this type of solutions in Figure 5.6, where the roadway merges 
at x = 0. For the case when the roadway diverges at x — 0, we can find similar 
solutions. 

From Figure 5.6, we obtain the boundary flux /q = f{U R ) for wave solutions of 
type 3. 

Type 4 When U R is in region OU^DO shown in Figure 5.2; i.e., 

f(U R ) < /([/*) = f(U L ), Pn/a R < p*/a* and a R < a*, (5.20) 

wave solutions to the Riemann problem are of type 4. These solutions consist of 
three basic waves with two intermediate states: U\ = (o,l, Pi\f(a L ,p 1 )=f{u 2 )) an d 
U2 = (or, P2\p 2 /a R =a)- Of these three waves, the left one (Ul, Ui) is a shock wave 
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with negative speed a = ^ 1 _ L ' < 0, the middle one (C/i,C/ 2 ) is a standing 
wave with zero speed, and the right one (U2,Ur) is a rarefaction wave with 
characteristic velocity Ai(o,p) > 0. 

From Figure 5.2, we can see that this type of solutions are admitted only when 
the roadway merges at x — 0. Here we present an example of this type of 
solutions in Figure 5.7. 

From Figure 5.7, we obtain the boundary flux /q = /(L^) for wave solutions of 
type 4. 

When Ul = (cll,Pl) is overcritical; i.e., Pl/ul > ct, where a is defined in (5.13), 
we denote the special critical point on 1-wave curve passing Ul as [/*; i.e., £/* = 
(oL)/0*|/9*/ai=a)- Thus, as shown in Figure 5.3, the [/-space is partitioned into three 
regions by three curves DU* = {a = a* = ai,0 < p < p*}, OU* = {0 < a < 
a*, p = aa} and C/*C = {a > o*, /(a, p) = /([/*). Related to different positions of the 
right state £/# in the t/-space, the Riemann problem for (5.5) with initial conditions 
(5.14) has the following six types of wave solutions. After discussion for each type of 
solutions we provide formula for calculating the associated boundary flux /q . 

Type 5 When Ur resides in region ABU^DA shown in Figure 5.3; i.e., 

f(U R ) < /(C/*), p R /a R < a and a R >a* = a L , (5.21) 

wave solutions to the Riemann problem are of type 5. These solutions con- 
sist of three basic waves with two intermediate states: U\ = [/* and Ui = 
(an, P2\f(u 2 )=f(u*))- Of these three waves, the left one (Ul,U±) is a rarefaction 
wave with negative characteristic wave velocity Ai(o, p), the middle one (U±, U2) 
is a standing wave and the right one (U2, Ur) is a rarefaction wave with positive 
characteristic velocity Ai(o, p). 

From Figure 5.3, we can see that this type of solutions are admitted only when 
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Figure 5.8: One example for wave solutions of type 5 for (5.5) with initial conditions (5.14) 






Figure 5.9: One example for wave solutions of type 6 for (5.5) with initial conditions (5.14) 
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Figure 5.10: One example for wave solutions of type 7 for (5.5) with initial conditions (5.14) 
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Figure 5.11: One example for wave solutions of type 8 for (5.5) with initial conditions (5.14) 
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Figure 5.12: One example for wave solutions of type 9 for (5.5) with initial conditions (5.14) 




Figure 5.13: One example for wave solutions of type 10 for (5.5) with initial conditions (5.14) 
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the roadway diverges at x = 0; i.e., clr > cll- Here we present an example of this 
type of solutions in Figure 5.8. 

From Figure 5.8, we obtain the boundary flux /q = f(Uo) for wave solutions of 
type 5. 

Type 6 When Ur resides in region BU^CBAs shown in Figure 5.3; i.e., 

f(U R ) > /(£/*), (5.22) 

solutions to the Riemann problem are of type 6. These solutions consist of three 
basic waves with two intermediate states: U\ = [/* and Ui = (cir, P2\f(u 2 )=f(u*))- 
Of these three waves, the left one (Ul, Ui) is a rarefaction wave with negative 
characteristic velocity Ai(o, p), the middle one (C/i, C/2) is a standing wave and 
the right one (U2, Ur) is a shock wave with positive speed a = ( R >ji( ±L . 

From Figure 5.3, we can see that this type of solutions may be admitted when 
the downstream traffic condition is undercritical or overcritical; However, they 
are admitted only when the roadway diverges at x — 0. Here we present an 
example of this type of solutions in Figure 5.9, where the downstream traffic 
condition is overcritical. For the case when the downstream traffic condition is 
undercritical, we can find similar solutions. 

From Figure 5.9, we obtain the boundary flux /q = /(t/2) for this type of wave 
solutions. Here we have the same formula as that for wave solutions of type 5. 

Type 7 When Ur resides in region CU^FUlEC shown in Figure 5.3; i.e., 

f(U L ) < f(U R ) < /([/,) and pr/clr > a, (5.23) 

wave solutions to the Riemann problem are of type 7. These solutions consist 
of two basic waves with an intermediate state U\ = (ox,, Pi\f(Ui)=f(u R ))- Of these 
two waves, the left one (Ul, Ui) is a rarefaction with negative characteristic 
velocity Ai(o, p), and the right one (U±, Ur) is a standing wave. 
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From Figure 5.3, we can see that the Riemann problem may admit this type 
of solutions when the roadway merges or diverges at x — 0. Here we present an 
example of this type of solutions in Figure 5.10, where the roadway diverges at 
x = 0; i.e., cir > aj_,. For the case when the roadway merges, we can find similar 
solutions. 

From Figure 5.10, we obtain the boundary flux /q = J'{Ur) for wave solutions 
of type 7. 

Type 8 When Ur locates in region FU^EOF shown in Figure 5.3; i.e., 

f(U R ) < f(U L ) < /([/*) and p R /a R > a, (5.24) 

wave solutions to the Riemann problem are of type 8. These solutions consist 
of two basic waves with an intermediate state U\ = (o>l, Pi\f(Ui)=f{u R ))- Of these 



two waves, the left on (Ul,Ui) is a shock with negative speed a 



f(U L )-f(Ui) 

and the right one (U±, Ur) is a standing wave. 

From Figure 5.3, we can see that the Riemann problem may admit this type 
of solutions when the roadway merges or diverges at x — 0. Here we present an 
example of this type of solutions in Figure 5.11, where the roadway diverges at 
x = 0; i.e., clr > a^. For the case when the roadway merges, we can find similar 
solutions. 

From Figure 5.11, we obtain the boundary flux /q = /(Ur) for wave solutions 
of type 8. Here we have the same formula as that for wave solutions of type 7. 

Type 9 When Ur resides in region DU^FGD shown in Figure 5.3; i.e., 

f(U L ) < f(U R ) < /(C/*), p R /a R < a and a R < a* = a L , (5.25) 

wave solutions to the Riemann problem are of type 9. These solutions consist of 
three basic waves with two intermediate states: U\ = (a L , Pi\f(Ui)=f(u 2 )) an d U 2 = 
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( a R, P2\p 2 /a R =a)- Of these three waves, the left on (Ul, Lq) is a rarefaction with 
negative characteristic velocity Ai(o, p), the middle one (Lq,C/ 2 ) is a standing 
wave, and the right one (U 2 , Ur) is a rarefaction with positive speed Ai(o, p). 

From Figure 5.3, we can see that this type of solutions are admitted only when 
the roadway merges at x = 0; i.e., a R < a^. Here we present an example of this 
type of solutions in Figure 5.12. 

From Figure 5.12, we obtain the boundary flux /q = f(U 2 ) for wave solutions 
of type 9. 

Type 10 When Ur resides in region GFOG shown in Figure 5.3; i.e., 

f(U R ) < f(U L ) < /(£/*)> Pn/a R < a and a R < a* = a L , (5.26) 

wave solutions to the Riemann problem are of type 10. These solutions consist 
of three basic waves with two intermediate states: Lq = (cll, Pi\f(Ui)=f(u 2 )) an d 
U 2 = ( a R, P2\ P2 /a R =a)- Of these three waves, the left one (Ul, U±) is a shock with 
negative speed, the middle one (Ui,U 2 ) is a standing wave, and the right one 
(U 2 , Ur) is a rarefaction wave with positive characteristic velocity \\(a, p). 

From Figure 5.3, we can see that this type of solutions are admitted only when 
the roadway merges at x = 0; i.e., clr < aL- Here we present an example of this 
type of solutions in Figure 5.13. 

From Figure 5.13, we obtain the boundary flux /q = f(U 2 ) for wave solutions 
of type 10. Here we have the same formula as that for wave solutions of type 9. 

5.3.2 Summary 

In the above subsection, we have discussed 10 types of wave solutions. For each type 
of solutions, the boundary flux /q is equal to one of the following four quantities: 
the upstream flow rate /(Ul), the downstream flow rate f(U R ), the capacity of the 
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No. 


left state Ul 


right state U R 


boundary flux /q 


1 


undercritical 


f(U R ) < f(U L ), a R > a*, p R /a R < a 


Wl) 


2 


undercritical 


f(U R ) > f(U L ) 


f(u L ) 


3 


undercritical 


f(U R ) < f(U L ), p R /a R > a 


f(u R ) 


4 


undercritical 


f(U R ) < f(U L ), p R /a R <a, a R < a* 


fmax 
JR 


5 


overcritical 


f(U R ) < fL ax , a R > a L , p R /a R < a 


fmax 
JL 


6 


overcritical 


f(u R ) > fr x 


fmax 
JL 


7 


overcritical 


f(U L ) < f(U R ) < fr\ PR/a R > a 


Wr) 


8 


overcritical 


f(U R ) < f(U L ), pn/a R > a 


Wr) 


9 


overcritical 


f(U L ) < f(U R ) < f? ax , p R /a R <a,a R < a L 


fmax 
JR 


10 


overcritical 


f(U R ) < f(U L ), p R /a R <a, a R < a L 


fmax 
JR 



Table 5.1: Solutions of the Boundary Fluxes 



upstream roadway /™ a:E and the capacity of the downstream roadway f R ax . For wave 
solutions of type 1 and 2, the boundary flux is equal to the upstream traffic flow rate; 
i.e., /q = /(Ul). For wave solutions of type 3, 7 and 8, the boundary flux is equal 
to the downstream traffic flow rate; i.e., /q = f(U R ). For wave solutions of type 4, 9 
and 10, the boundary flux is equal to the capacity of the downstream roadway; i.e., 
/o = fR ax - F° r w ave solutions of type 5 and 6, the boundary flux is equal to the 
capacity of the upstream roadway; i.e., /q = f™ ax . In Table 5.1, the boundary fluxes 
are listed for the 10 types of wave solutions to the Riemann problem, as well as the 
conditions when the Riemann problem admit those solutions. 

Note that when a^ = a R ; i.e., when (5.4) becomes a homogeneous LWR model, 
wave solutions and the solutions of the boundary fluxes provided here are the same 
as those for the homogeneous LWR model. 

The Riemann problem for an inhomogeneous LWR model was also studied by 
Lebacque (1995). He discussed the Riemann problem for (5.5) with general inho- 
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Conditions 


Solutions by Lebacque 


Types 


Our solutions 


cll < a R , U L uc, U R uc 


Wl) 


1 


Wl) 


a-L < a R , U L uc, U R oc 


mm{f(U L )J(U R )} 


2 or 3 


f(U L ) or f(U R ) 


«l < a R , U L oc, U R uc 


fmax 
JL 


5 or 6 


fmax 
JL 


a-L < cir, U L oc, U R oc 


mm{fr x J(U R )} 


6, 7 or 8 


f max Qr /([/r) 


ol > clr, U l uc, U r uc 


mm{f% a *,f(U L ) 


1 or 4 


w L ), fr x 


Ol > Ofl, U L UC, f/ij oc 


mm{f(U L )J(U R )} 


2 or 3 


f(U L ) or f(U R ) 


OL > CLR, U L OC, C/fl uc 


fmax 
JR 


9 or 10 


fmax 
JR 


ClL > ClR, U L OC, Ur OC 


Wr) 


7 or 8 


Wr) 



Table 5.2: Comparison with Lebacque's results 



mogeneity, and presented empirical solutions for it. He categorized the solutions 
according to two criteria. The first criterion is to consider the relationship between 
the upstream capacity and the downstream capacity. For the roadway with variable 
number of lanes, this criterion is equivalent to considering the relationship between 
the number of lanes of the upstream and downstream roadway, since the roadway 
with greater number of lanes has larger capacity. The second criterion is to con- 
sider whether the upstream and downstream traffic conditions are undercritical or 
overcritical. With these criteria, he discussed 8 types of waves solutions to the Rie- 
mann problem and obtained the formula for the boundary flux related to each type 
of solutions. The conditions for those types of wave solutions as well as the formulae 
related to those types of solutions are listed in Table 5.2, in which oc and uc stand for 
overcritical and undercritical respectively. Under each of those conditions, the Rie- 
mann problem may admit different types of solutions discussed the above subsection 
5.3.1. The types of solutions and our related formulae for the boundary flux are also 
presented in Table 5.2. From this table, we can see that our results are consistent 
with those provided by Lebacque, although the Riemann problem is solved through 



CHAPTER 5. The Inhomogeneous LWR Model and Its Numerical Solutions 103 

different approaches. 

The consistency of our results with existing results can also be shown by intro- 
ducing a simple formula for the boundary flux. If we define the upstream demand 

as 

{f(UA plI^l < ol 
(5.27) 
fT ax PL/a L > a 

and define the downstream supply as 

f R = ( ST * PB/a « < ° (5.28) 

[ f(U R ) pn/a R > a 

then the boundary flux can be simply computed as 

/* = mini/*,/*}. (5.29) 

Note that /£ = /(£/*). Formula (5.29) was also provided by Daganzo (1994, 1995) 
and Lebacque (1995). 

5.4 Numerical solution method 

Since the inhomogeneous LWR model can be written in a conservation form (5.5), 
Godunov's method is efficient for its numerical solutions. In this section, we describe 
the Godunov method for solving (5.5) with general initial and boundary conditions. 
In a Godunov's method for (5.5), the roadway is partitioned into N zones and a 
duration of time is discretized into M time steps. In a zone i, we approximate the 
continuous equation (5.4) with a finite difference equation 

At + A^ " U ' (5 ' dUj 

where p™ 1 denotes the average of p in zone i at time step m, similarly p™ +1 i s ^e 
average at time step m + 1; f* + i/ 2 denotes the flux through the upstream boundary 
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of zone i, and similarly f* +1 / 2 denotes the downstream boundary flux of zone %. In 
(5.30), the boundary flux f*_i/ 2 is related to solutions to a Riemann problem for (5.5) 
with the following initial conditions: 

{U™! X < Xi-i/2 
(5.31) 
U™ x> Xi- 1/2 

The wave solutions to the Riemann problem and the associated formula for the bound- 
ary flux have been discussed in section 5.3. Then according to (5.30) we can find p at 
time step m + 1, given traffic conditions at time step m and values of p at the road 
boundaries. 

5.5 Conclusions 

In this chapter, we find that the inhomogeneous LWR model (5.5) exhibits nonlinear 
resonance. For this system, we discuss the Riemann problem and find 10 types of 
wave solutions and the formula for the boundary flux related to each type of wave 
solutions. Then we obtain a general formula for the boundary flux with the defini- 
tion of upstream demand and downstream supply. We also conclude that the results 
obtained here are consistent with those in existing literature. Since the inhomoge- 
neous LWR model can be written in a conservation form (5.5), Godunov's method is 
efficient for its numerical solutions. Godunov's method for the inhomogeneous LWR 
model is briefly described in section 5.4. 

In this chapter, we consider the inhomogeneity as variable number of lanes. Sim- 
ilarly we can solve the inhomogeneous LWR model with other inhomogeneities, by 
changing the functions v(a,x) and f(a,x) corresponding to the new inhomogeneities. 
For all kinds of inhomogeneity, equation (5.29) is always the general formula for 
computing the boundary flux. 

In this chapter, we consider an inhomogeneous LWR model, which is a first-order 
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traffic flow model. However, we want to point out that the method provided here can 
be extended to solve inhomogeneous higher-order traffic flow models such as Zhang's 
model discussed in Chapter 3 and the PW model discussed in Chapter 4. 



Chapter 6 

A First-Order Multi-Commodity 
Model and Its Numerical 
Simulations 



6.1 Introduction 

All the models discussed in previous chapters are for link flows, in which all vehicles 
have the same contribution to the flow dynamics. When one considers a network flow, 
he/she has to deal with vehicles of different origins, destinations, or other attributes. 
These attributes of vehicles play a role in determining their choice of routes, and 
therefore affect the flow dynamics. Hence, to model a network flow, one has to take 
these attributes into account and disaggregate traffic flow into different components. 
Such a network flow model considering disaggregated traffic flow is called a multi- 
commodity model. 

In literature, there have been several models for multi-commodity flow, including 
the model by Vaughan, Hurdle and Hauer (1984), Jayakrishnan's model (1991) and 
Daganzo's cell transmission model (1994; 1995). 
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Vaughan, Hurdle and Hauer (1984) suggested a continuous model consisting of a 
'local equation' and a 'history equation'. In their model, the basic element of traffic 
flow is a vehicle, and for each vehicle a label for its trajectory is introduced. In this 
model, the local equation ensures traffic conservation, and the history equation com- 
putes the trajectory of each vehicle. It has been shown that this model is consistent 
with the traditional LWR model at the aggregate level. 

Jayakrishnan (1991) introduced another multi-commodity model. This is a dis- 
crete model, in which each link is partitioned into a number of zones. In each zone, 
vehicles close to each other and with the same origin, destination or other common 
commodity-specific characteristics are considered as a "macroparticle" . To determine 
the position of a macroparticle at next time step, Jayakrishnan considered its travel 
speed and the length of the zone where it stays. This model doesn't always preserve 
traffic conservation; i.e., this model may not be consistent with the LWR model, since 
the macroparticles are considered separately. 

Daganzo (1994; 1995) introduced a multi-commodity model based on his cell 
transmission model. In this discrete mult i- commodity model, traffic flow in a zone is 
also disaggregated into macroparticles. In every zone, the macroparticles are ordered 
according to the waiting time. In this order, the macroparticles are moved, and those 
macroparticles can proceed to a downstream zone if their waiting times are greater 
than a threshold minimum waiting time. In Daganzo's model, the minimum waiting 
time is set to the waiting time of a macroparticle if the total number of vehicles 
which are in the same zone and in front of this macroparticle is equal to the number 
of vehicles moving into the downstream zones. In this model, the number of vehicles 
moving into the downstream zones are computed as the boundary flux times the 
length of a time step based on solutions to the Riemann problem for the LWR model, 
which has been discussed in Chapter 5. 

Of these three models, Jayakrishnan's and Daganzo's are simpler in computation 
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since a macroparticle instead of a vehicle is considered as the basic element of traffic 
flow. 

Both the model by Vaughan et al. and Daganzo's model consider the dynamics 
of a network flow at two levels: aggregate level and disaggregated level. At aggregate 
level, vehicles in a network flow is considered to be identical particles. At this level, 
both models are consistent with the LWR model, and therefore they satisfy traffic 
conservation. In the model by Vaughan et al., the continuous LWR model is used, 
and while Daganzo's model uses the discrete form of the LWR model. 

In these three models the vehicles are ordered. In the model by Vaughan et 
al. vehicles are ordered according to the trajectory labels, Jayakrishnan ordered 
macroparticles according to location and Daganzo ordered the macroparticles accord- 
ing to time. They all assume that the vehicles observe the First-In-First-Out (FIFO) 
discipline when they travel on the roadway. In particular, the model by Vaughan 
et al. assumes the vehicle trajectories don't cross each other at any time t and any 
location x, Jayakrishnan's model assumes that vehicles always keep the order in lo- 
cation, and Daganzo's model assumes that vehicles always keep the order in time. Of 
all these models, the one by Vaughan et al. uses the full information of order in time 
and location, while the other two use part of the order information of either time or 
location. 

In this chapter we introduce a new multi-commodity model. Our model is based 
on the concept of "macroparticle", and has also a two-level structure. In our model, 
we use the discrete form of the LWR model for traffic flow at the aggregate level. In 
our new model, we will interpret the FIFO discipline in a new way so that we can 
simplify the dynamics at disaggregated level. 

This chapter is organized as follows. In section 2, we discuss the new multi- 
commodity model and define the network structure, data structure and program 
structure. In section 3, we present some numerical simulations. In section 4, we 
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discuss the future plans on this topic. 

6.2 A two-level multi-commodity model 

The dynamics of a network flow is considered at two levels: aggregate level and 
disaggregated level. 

For traffic flow at the aggregate level, we use the discrete form of an inhomoge- 
neous LWR model. As discussed in Chapter 5, the inhomogeneous LWR model can 
be written as 

(H + f(a,p) x = 0, (6.1) 

in which a(x) is an inhomogeneity factor. By partitioning each link into N zones, and 

discretizing the time interval into M time steps, we obtain the Godunov-type finite 

difference equation for (6.1): 

j+i j tj* _ fj* 

At Ax y ' 

where Ax is the length of zone i, At is the time from time step j to time step j + 1, 
and the choice of -^ is governed by the CFL condition. In equation (6.2), pi is the 
average of traffic density p in zone i at time step j , similarly pi is the average of p 



at time step j + 1; fl^ x i 2 is the flux through the upstream boundary of zone % from 



c3* 

/2 1= 111C UUA WllUUfcll 

ej* 



time step j to time step j + 1, and similarly f i+l / 2 is the downstream boundary flux. 
Given traffic conditions at time step j, we can calculate the traffic density in zone % 
at time step j + 1 as 

Pf 1 = P? + ^(/£i/2-/£i/ 2 ), (6-3) 

where computation of the boundary fluxes fWi 2 an d //+i/ 2 nas been discussed in 
Chapter 5. 

Defining N- = p > i Ax as the number of vehicles in zone i at time step j, Nf + = 
pl + Ax as the number of vehicles at time step j + 1, F^_ 1/2 = Atf{p 3 *_ l/2 ) as the 
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number of vehicles flowing into zone i from time step j to j + 1, and F? +1 , 2 as the 
number of vehicles flowing out of zone i, equation (6.3) can be written as: 

N j+1 = N j + F j , - F j , (6 A) 

ly i ly i ~ 1 i— 1/2 1 1+1/2' V U, V 

which is a conservation form. 

For all the link flows in a network, equation (6.4) is an efficient model at aggregate 
level, which captures the dynamics and preserves the conservation of traffic flow. 
However, how to model the dynamics for the flows at merges and diverges is still an 
open question. 

There have been two empirical treatments of flows at merges by Daganzo and 
Lebacque. For a merge with K upstream zones, Daganzo (1995a) solve the flux 
through the merge using traffic supply of the downstream zone, and the sum of traf- 
fic demands of these K upstream zones as traffic demand. This flux is distributed 
to all the upstream zones according to some pre-defined distribution fractions. This 
treatment has an underlying assumption that the upstream flows have the same con- 
tributions to the merge. This assumption limits the application of this treatment for 
a merge where the upstream zones consist of mainline highway and on-ramps, since 
on-ramps have lower priority than mainline highway. Another approach suggested by 
Lebacque (1995) first splits traffic supply of the downstream zone into K parts ac- 
cording to some pre-defined fractions, and forms a pair of demand and supply for each 
of the K upstream zones. Then for each pair of demand and supply, a flux through 
the merge can be solved. The sum of these K fluxes, considered by Lebacque, is then 
the total flux through the merge. In our model, we will use Lebacque's method to 
treat the merges. 

For traffic flow at diverges, the approaches suggested by Daganzo and Lebacque 
can still be used. In Daganzo's method, the traffic supply for a diverge with K 
downstream zones is the sum of K traffic supplies. Then the flux through the diverge 
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is computed with this traffic supply and the traffic demand of the upstream zone and 
distributed to the downstream zones according to pre-defined fractions. In Lebacque's 
method, for each of K downstream zones, the traffic demand is a fraction of the 
upstream traffic demand. With this traffic demand and its traffic supply, the flux is 
computed. The sum of the K fluxes is the total flux through the diverge. In our 
model, we use a treatment which is adjusted from Lebacque's. Here, we split the 
upstream demand based on the route choices instead of some pre-defined fractions; 
i.e., the demand for each downstream zone is equal to the number of vehicles choosing 
that zone. This treatment may better model the route choices at diverges. 

At disaggregated level, we partition vehicles in a zone into macroparticles. In 
a macroparticle, vehicles are close to each other and have the same disaggregated 
information such as destination. The macroparticles in a zone are ordered in location 
as a queue, with those at downstream as the head of the queue and those at upstream 
as the tail of the queue. Since the CFL condition guarantees that no macroparticle 
can cross a zone in one time step, our question is whether a macroparticle can be 
moved into a downstream zone and what position it is in that zone. 

To determine the movement of macroparticles, we introduce a new concept - " 
boundary connector". A boundary connector stands for a boundary between linked 
zones in a link of a traffic network. It is not a physical entity, doesn't have length 
and cannot store vehicles. With this concept, zones are not considered to be linked 
to each other directly, but are linked to boundary connectors. Therefore, in a traf- 
fic network, zones reflect traffic conditions, and boundary connectors determine the 
network structure. 

For each boundary connector, we determine the movement of macroparticles at 
every time step as follows. From each upstream zone connected to the boundary con- 
nector, we pick out macroparticles starting from the head of the queue of macropar- 
ticles and store them in the boundary connector temporarily, until the total number 
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of vehicles moved is equal to the flux from that zone, which is determined by the 
computation at aggregate level. When there are more than one zones connected to 
the boundary connector, the macroparticles moved into the connector are ordered 
according to the times they enter. However, in reality, we do not know the exact 
order of their arrival times and a random process may be used to order them. Then 
every macroparticle is moved into the downstream zone leading to its destination. 
If a macroparticle has more than one downstream zones to choose, corresponding 
fractions may be used. After a macroparticle is moved into a zone, it is attached to 
the tail of the queue in that zone. In our model, the FIFO principle is also applied. 
However, this principle is interpreted in a different way than in Vaughan et al., Da- 
ganzo or Jayakrishnan. Here, FIFO means those macroparticles in front (location) 
in a zone are ahead (time) in a boundary connector, and similarly those ahead in a 
boundary connector are in front in a zone. 

The model suggested here has been implemented for a simple network. In the 
remaining part of this chapter, we design the network, data and program structures 
for it and carry out some numerical simulations. 

6.3 Network, data and program structures for a specific traf- 
fic network 

In this section a simple traffic network, shown in Figure 6.1, is considered. This 
network consists of a mainline highway, one on-ramp and one off-ramp. Therefore, 
there are two origins and two destinations in this network. For modeling purpose, the 
mainline highway is partitioned into 20 zones, labeled from 1 to 20, and 21 boundary 
connectors are used to connect these 20 zones, highway origin, highway destination 
and two ramps. Of these 21 boundary connectors, boundary connector 1 is the 
upstream highway boundary and boundary connector 21 is the downstream highway 
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boundary. Besides 20 highway zones, we use artificial zone to denote highway origin 
1, zone 21 to denote highway destination 1, zone 22 to denote on-ramp, and zone 23 
to denote off-ramp. 

The data used to reflect network traffic flow are shown in Figure 6.2. At aggre- 
gate level the data include: 

1. Traffic measurements of zone i at time step j - p?, v\. The number of vehicles 
N? = p^Axi, where Axi is the length of zone i. These measurements are given 
initially and are updated at each time step. The traffic conditions of the origin 
zones 1 and 22, the destination zones 21 and 23 may be given in certain type of 
boundary conditions. 

2. Fluxes through boundary i + 1/2 between time step j and time step j + 1. We 
solve the Riemann problem to get p J -l 1 / 2 , vJ li/o- The number of vehicles through 



the boundary F? +1 , 2 = p 3 j ^ 1 , 2 v 3 i ^ 1 , 2 AP , where At 3 is length of the time interval 



t+l/2' t+1/2' 
i+1/2 ~~ l J i+l/2 u i+l/7 

from j to j + 1. The fluxes at boundary 1 and boundary 21 may be given with 
certain type of boundary conditions. 

At disaggregated level the data include 

1. Array of queues. Each zone is a dynamically allocated queue. Each node of a 
queue stands for a macroparticle. The parameters for a macroparticle include its 
destination, the number of vehicles in it and a pointer to its upstream macropar- 
ticle. The parameters of a zone include the number of macroparticles in this 
zone, the number of vehicles to each destination, the total number of vehicles in 
the zone and some other related information of the zone. 

2. Array of boundary connectors. Each boundary connector stores the upstream 
and downstream zones that are connected to it. For different types of boundary 
connectors, e.g., merge or diverge, the treatments are different, therefore the 
type of a boundary connector is a parameter. Another parameter is the number 
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Figure 6.2: Data structure of traffic network 
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of vehicles passing through a boundary connector at a time step. Under some 
conditions, these fluxes through some specific boundaries may be set, for exam- 
ple, at boundary connector 2, the flux from zone 22 may be set to a value if this 
on-ramp is metered at a rate. 

The program flow chart for computing aggregate flows is shown in Figure 6.3. 
This program consists of the following operators: 

1. The Riemann solver. The Riemann solver is used to calculate the fluxes through 
the boundaries. The solver can be of first- or second-order. 

2. Updating of traffic conditions. This is done according to the finite difference 
equations. For the LWR model, traffic conditions are updated according to 
(6.4). 

When disaggregated traffic flow is considered, the program flow chart is shown 
in Figure 6.4. This program consists of the following operators: 

1. Initialization of traffic condition, which provides the initial traffic measurements 
to all the zones and macroparticles. One simple way to start is to assume that 
the network has no traffic in the beginning. 

2. Initialization of network structure, which provides values for the parameters of 
each boundary connector. 

3. The Riemann solver. The solver is used to compute the aggregate traffic flux. 
It gives the number of vehicles through the boundary connector during a time 
interval. 

4. Updating of traffic conditions around a boundary connector. At each time step, 
the traffic conditions of those zones connected to a boundary connector are up- 
dated with the flux through the boundary. We retrieve macroparticles from 
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upstream zones, and distribute them to downstream zones according to the des- 
tination information carried by each macroparticle. This process stops when 
the number of vehicles passing the boundary is equal to the flux calculated by 
the Riemann solver. The retrieving from upstream zones is done by the func- 
tion detach; and the distributing to downstream zones is done by the function 
attach. The function communicate governs these two operators according to 
the O-D information carried by macroparticles. This function uses certain merg- 
ing rules and route choice rules which are used at aggregate level. Some ad-hoc 
treatment on these rules was discussed in previous section. 

5. Providing demands at origin zones. The demand at each origin zone is given as 
a queue of well-ordered macroparticles, although in reality this information is 
usually unknown. 

6. Computing output at destination zones. The traffic conditions of each destina- 
tion zone are the outputs of interest. 

6.4 Numerical Simulations 

In this section, we carry out one numerical simulation for the network shown in 
Figure 6.1, in which zone 2 and zone 18 have 4 lanes, zone 22 and zone 23 have 1 
lane; all the other zones have 3 lanes. 

The LWR model is used for the aggregate flow with Newell's fundamental diagram 
f(p) = P v f[L — exp( — (1 — pj/p))}, shown in Figure 6.5. The free flow speed Vf = 60 
mph, the wave speed of jam density Cj = —10 mph and the jam density pj = 250 
vpm. The characteristic speed of the LWR model A* — Vf * (1 — e [Cj ^ v ^ 1 ~ p ^ p ^) — 
CjPj/ pe^ Cj ^ Vf<yl ~pjl p \ The length of each zone Ax = 0.6 mile; the length of a time step 
At = 30 sec. The CFL number is 60 • 5/600/0.6 = 0.8444. The maximum number 
of vehicles in a zone is 250x0.6=150. The maximum number of vehicles through a 
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Figure 6.5: Newell's Fundamental Diagram 
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boundary in At is about 12 vehicles, when p = 60vpm, / = 1476vph. 

We assume origin zone is jammed with a repeating sequence of platoons with 
25 vehicles to destination 1 followed by 10 vehicles to destination 2. Origin zone 22 
is assumed to be jammed with a repeating sequence of platoons with 5 vehicles to 
destination 1 followed by 2 vehicles to destination 2. We assume destination zone 21 
has the same number of vehicles as zone 20, and destination zone 23 has the number 
of vehicles to destination 2 in zone 18. 

Using the merging rule introduced by Lebacque (1995), two Riemann problems 
are solved at boundary connector 2 and the boundary flux is calculated. At boundary 
connector 19, two fluxes to destination zone 23 and zone 19 are obtained by solving 
two Riemann problems. The sum of them is the total flux through the boundary. 

We get the following numerical results: 

1. The output of two destinations are given in Figure 6.6. 

2. The fluxes through all the boundary connectors are given in Figure 6.7 

3. The number of vehicles of all the zones are given in Figure 6.8 

4. The densities and velocities in zones 1 to 21 are given in Figure 6.9 

According to the results, a rarefaction wave still forms. However, when traffic 
flow moves to boundary connector 19, all the vehicles to destination 2 move out of 
the mainline of highway. (Refer to the below part of Figure 6.8) Zone 2 becomes 
jammed after about 60 x 30 seconds, since traffic has to merge from 4 lanes in zone 
2 to 3 lanes in zone 3. These are consistent to experiences. However we do not have 
field observation data to justify our model at this time. 
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6.5 Conclusions 

In this chapter, a new multi-commodity model is developed. In this model, by inter- 
preting the FIFO principle in a new way and introducing the concept of boundary 
connectors, we develop a more efficient approach for disaggregated flows. Therefore, 
this new multi-commodity model is promising for modeling complex network flows. 

For traffic flows at aggregate level, higher-order traffic flow models may be com- 
bined in order to capture more dynamics. For example, the PW model or Zhang's 
model discussed in previous chapters can be combined for link flows. 

However, how to model the dynamics at merges or diverges is still an open prob- 
lem, which is much harder than to model link flows. 

Since this model is discrete, and network, data and program structures have been 
discussed in detail, it is ready for practical test. For example, this multi-commodity 
model can be used in simulating network traffic conditions in order to verify it validity. 

Some possible applications of this network flow model can be seen now. One 
application is for developing better ramp metering methods. Another application is 
to determine dynamic assignment in a complex traffic network. 



Chapter 7 



CONCLUSIONS 



This chapter provides concluding remarks on the research effort reported in this thesis. 
The chapter starts with an overall conclusion in section 7.1; In section 7.2 we discuss 
our contribution to traffic flow models and their numerical solutions. In the last 
section we discuss the future research directions in some areas related to this topic. 

7.1 Overall Conclusions 

In this research we studied five traffic flow models: the LWR model, Zhang's model, 
the PW model, the inhomogeneous LWR model and the mult i- commodity model. For 
these models, we discussed their analytical solutions, developed numerical solution 
methods and carried out numerical simulations. 

All these models preserve conservation of traffic. Since the homogeneous LWR 
model is the simplest model preserving traffic conservation, it is the basis to under- 
stand the wave solutions of these models and develop numerical solution methods 
for them. The LWR model is a first-order hyperbolic conservation law'. Shock and 
rarefaction waves are the basic solutions to a Riemann problem for such a conser- 
vation law. From the solutions to the Riemann problem, we can easily compute the 
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boundary averages of traffic density and boundary fluxes, which are used in the nu- 
merical solution methods - Godunov's methods. We solve the homogeneous LWR 
model with a first-order Godunov method and our numerical results are consistent 
with theoretical predictions. 

In chapter 3 we discussed a non-equilibrium traffic flow model - Zhang's model, 
which is a second-order hyperbolic system of conservation laws with a source term. 
Due to the difficulty in solving the Riemann problem for a system with a source term, 
we studied wave solutions to the Riemann problem for the homogeneous version of 
this model without considering the source term. The wave solutions are much more 
complicated than those to the Riemann problem for the LWR model. From these wave 
solutions, we computed the boundary averages of p and v, and developed a first-order 
and second-order Godunov methods for Zhang's model. The performance of these 
methods are examined in Section 3.4. We found that the second-order Godunov 
method performs better than a first-order method, however, its converge rate is 1, 
instead of 2 as expected and this is believed to be caused by the effect of the source 
term. 

In Chapter 4 we discussed another non-equilibrium second-order traffic flow model 
- the PW model. We discussed wave solutions to the Riemann problem for the 
homogeneous version of the PW model, as well as wave solutions to the Cauchy 
problem for the PW model with a source term. For the Cauchy problem, we found 
that the characteristics of a I-rarefaction wave are approximated by parabolic curve. 
However, numerical results didn't show significant improvement in solutions to the 
PW model by solving Cauchy problems than by solving Riemann problems. For 
the PW model, we studied a first-order Godunov method, a second-order Godunov 
method, Pember's method, fractional method and LeVeque's method. All the higher- 
order methods don't show significantly better convergence rates than the first-order 
method, due to the effect of the source term. In LeVeque's method, we interpolated 
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the given p and v with solutions of the standing wave in each cell. However, this 
method was proved to under-estimate the effect of the source term and is not suitable 
for the PW model when solutions are far from equilibrium states. Different from 
Zhang's model, the PW model is unstable in certain density/speed regions. Therefore 
we examined the stability of the PW model with a first-order Godunov method in 
Section 4.4.1. The results obtained are consistent with theoretical predictions. 

In Chapter 5 we studied the inhomogeneous LWR model. By introducing a profile 
of inhomogeneities the inhomogeneous LWR model can be written as a 2 x 2 non- 
strictly hyperbolic system of conservation laws. Then in Section 5.3 we rigorously 
solved the Riemann problem for the inhomogeneous LWR model for a roadway with 
variable number of lanes. We found that the results are consistent with those empirical 
results found in literature. However our method is easier to be extended to higher- 
order inhomogeneous models. 

The multi-commodity model discussed in Chapter 6 is different from other mod- 
els since it has a two-level structure and uses a discrete form of the LWR model 
at aggregate level. In Chapter 6, we made clear of the two-level structure of all 
multi-commodity models. By introducing boundary connectors, we suggested a more 
efficient method for computing the dynamics at disaggregated level. We also designed 
the network, data and program structures for a specific network in Section 6.3, and 
in Section 6.4, we presented some numerical simulation results which are consistent 
with theoretical expectations. 

Through discussions on different types of models and extensive numerical simula- 
tions, this thesis builds a solid base for modeling traffic flow macroscopically. It helps 
to understand the wave properties of the continuum models. The numerical methods 
discussed in this work will be applicable to all applications related with those models. 
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7.2 Research Contributions 

In this subsection, we discuss about our contributions in this thesis. 

Godunov methods have been well developed for studying hyperbolic systems of 
conservation laws. However they have not been extensively applied to the study 
of traffic flow problems. In this research, we made many efforts to apply Godunov 
methods and their variants to solve different traffic flow models. Our research shows 
that, to the continuum traffic flow models, Godunov methods are as useful as to other 
fluid dynamics. 

In traditional Godunov methods, boundary averages are computed from solutions 
to a Riemann problem for a homogeneous system. In this research (Section 4.2.2), 
we tried to compute boundary averages from solutions to a Cauchy problem for the 
system with a source term. Although this new approach in computing boundary 
averages doesn't appear to improve the numerical accuracy significantly, it may help 
people to design Godunov-type methods which are more suitable for systems with 
source terms. 

The inhomogeneous LWR model was considered as a 2 x 2 nonlinear resonant 
system in this thesis (in chapter 5), and rigorous solutions to the Riemann problem 
are studied based on existing theories, which were not known to the area of traffic 
flow models. This new approach helps us in understanding the wave phenomenon on 
an inhomogeneous roadway, and also brings us a new choice in solving inhomogeneous 
traffic flow models. 

In chapter 6, we made the two-level structure of mult i- commodity models clearer, 
and presented a way to deal with traffic flows at different levels more efficiently. The 
discussion on the network, data and program structures for traffic networks in this 
research is another contribution. 
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7.3 Future Research 

In this research, we made good progress in understanding some of the well-known 
continuum models. However, many topics discussed here still require further study. 
For example, an inhomogeneous model can be introduced for non-equilibrium flow 
and for multi-commodity flow, whose wave properties could be different than those 
in the homogeneous models. 

Modeling of traffic flow is becoming increasingly important. However, develop- 
ment of new realistic traffic flow models is still a challenging task. For example, there 
are few models that accurately capture the dynamics at merges and diverges or the 
interaction between different lanes on a multi-lane highway, a problem briefly touched 
in Chapter 6. 

Besides developing new models, another major challenge is the validation of traffic 
flow models that have been proposed, such as the PW model and Zhang's model. It's 
not simply a matter of computing the results and comparing with observed data, a 
deeper understanding of the wave phenomena found in traffic flows and how to model 
them are critical to this endeavor and should be further pursued. 

Although Godunov type of methods have been proved to be useful in solving 
traffic flow model, more efficient Godunov methods are needed to handle the diverse 
characteristics of these models. Another direction of research in computation may 
come from the development of parallel algorithms to simulate traffic flow in large 
networks. 
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